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1.  INTRODUCTION 


1.1  Summary 

It  is  known  that  large  space  structures  will  be  subjected  to 
thermomechanical  loadings  and  environmental  conditions  which  are  likely  to 
degrade  the  constitutive  properties  of  the  structural  materials,  thus  leading 
to  possible  failure  of  these  vehicles.  Therefore,  it  is  desirable  to  develop 
new  analytical  models  which  are  capable  of  accounting  for  these  degraded 
properties  so  that  design  procedures  can  be  improved.  There  are  three 
important  aspects  of  such  an  effort:  selection  and  development  of 
constitutive  models  applicable  to  large  space  structures,  construction  of 
analytical  models  and  experimentation  to  determine  the  precise  nature  of  the 
material  parameters  to  be  utilized  in  the  analytical  model.  These  three 
components  of  the  research  must  be  tied  together  into  a  single  concise  package 
in  order  to  obtain  a  useful  model. 

This  research  project  is  a  three  year  effort  to  develop  an  analytic  model 
capable  of  predicting  the  response  of  space  structures  with  degrading  material 
properties  under  quasi-static  as  well  as  dynamic  cyclic  thermomechanical 
loading  conditions.  The  research  was  funded  by  the  Air  Force  Office  of 
Scientific  Research  under  contract  no.  F49620-83-C-0067. 

1.2  Statement  of  Work 

Models  have  been  developed  for  predicting  the  thermomechanical  response 
of  large  space  structures  to  cyclic  transient  temperature  loading 
conditions.  The  research  was  conducted  in  the  following  stages: 

1)  selection  and  specialization  of  thermomechanical  constitutive 
equations  to  be  utilized  in  the  analysis  of  large  space  structures; 


2)  construction  (where  necessary)  of  coupled  energy  balance  equations 
(modified  Fourier  heat  conduction  equations)  applicable  to  the 
constitutive  models  selected  in  item  1); 

3)  casting  (where  necessary)  the  resulting  field  laws  into  coupled  and 
uncoupled  variational  principles  suitable  for  use  with  the  finite 
element  method; 

4)  finite  element  discretization  of  the  variational  principles  for  large 
space  structures; 

5)  experimentation  to  determine  material  properties  to  be  utilized  in  the 
constitutive  models;  and 

6)  parametric  studies  of  the  quasi -static  and  dynamic  response  of  large 
space  structures  undergoing  thermomechanical ly  and  environmentally 
degraded  material  properties. 


The  experimental  effort  (discussed  in  5)  was  supported  in  part  by  DOD 
equipment  grant  no.  841542.  The  total  research  effort  outlined  above  spanned 
a  period  of  three  years.  The  following  section  details  results  obtained 


2.  RESEARCH  COMPLETED  TO  DATE 


2. 1  Summary  of  Completed  Research 

The  following  tasks  have  been  completed  during  the  contract  period: 


1)  literature  survey; 

2)  selection  of  constitutive  equations  for  thermoviscoplastic  metals  at 
elevated  temperatures  and  polymeric  composites  with  thermomechanical 
load  induced  damage  at  temperatures  below  the  glass  transition 
temperature; 

3)  construction  of  a  coupled  energy  balance  equation  for 
thermoviscoplastic  metals; 

4)  casting  of  field  laws  for  the  material  discussed  in  2)  into  a  one¬ 
dimensional  finite  element  computer  code  with  two-way  thermomechanical 
coupling; 

5)  parametric  studies  using  the  model  developed  in  4)  to  determine  the 
thermomechanical  response  of  representative  metallic  space  structures 
with  degraded  material  properties; 

6)  development  of  generalized  constitutive  equations  for  metal  matrix 
composites  with  distributed  damage; 

7)  experimentation  to  determine  material  parameters  for  the  model 
developed  in  item  6) ; 

8)  development  of  algorithms  for  composite  truss-like  space  structures 
with  damage  induced  and  spacially  variable  stiffness  loss; 

9)  parametric  studies  for  graphite/epoxy  composite  space  structures  using 
’  item  8); 

10)  development  of  bounding  techniques  for  hysteretical ly  induced 
temperature  rise  in  thermoviscoplastic  space  structures; 

11)  development  of  an  analytic  method  for  modeling  beam-like  structural 
components  with  damage  Induced  stiffness  loss;  and 

12)  development  of  a  finite  element  model  for  composite  beam-like  space 
structures  with  elastic  material  properties  and  subjected  to  solar 
flux  heating  and  radiation  boundary  conditions. 


2.2  Literature  Survey 

A  detailed  literature  survey  has  been  completed  as  part  of  the  research 
effort  [1].  This  report,  entitled  "Large  Space  Structures  Technology:  A 
Literature  Survey,"  was  included  in  the  first  annual  technical  report. 
Briefly,  the  report  details  recent  advances  in  the  areas  of  materials, 
structural  solution  techniques,  damping,  and  preliminary  design  and 
experiment.  The  results  of  this  survey  indicate  that  very  little  research  is 
available  on  the  effect  of  material  property  degradation  on  large  structural 
response. 


2.3  Selection  of  Constitutive  Equations 

Candidate  material  models  have  been  selected  for  metals  at  elevated 
temperatures  and  polymeric  composites  below  the  glass  transition 
temperature.  These  are  detailed  below. 


2.3.1  Metals  at  elevated  temperatures  are  currently  modeled  using  continuum 
mechanics  with  internal  state  variables  (ISV's)  [1-5],  wherein  the  stress- 
strain  relation  is  of  the  form  (for  inf initesminal  strains): 
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where  is  the  stress  tensor,  ekl  is  the  strain  tensor,  D .  jki 
elastic  modulus  tensor,  is  the  inelastic  strain  tensor. 
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thermal  strain  tensor.  In  addition. 
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where  f.j  and  Q^.  are  appropriate  functions  of  state,  T  is  the  temperature, 
and  are  a  set  of  second  order  tensor  valued  internal  state  variables 

modeling  dislocation  arrangement  dislocation  density,  intergranular  damage. 


Although  it  has  been  demonstrated  that  numerous  models  fall  within  the 
above  framework  [61,  the  special  cases  of  equations  (1)  through  (3)  utilized 
thus  far  are  a  classical  plasticity  model  developed  by  Allen  and  Haisler  17,81 
(see  Appendix  6.2),  a  single  internal  state  variable  viscoplastic  model 
developed  by  Cernocky  and  Krempl  [9,10]  (see  Appendix  6.1),  and  a  two  internal 
state  variable  viscoplastic  model  developed  by  Bodner,  et  al.  [11-141  (see 
Appendices  6.1,  6.3,  6.4,  and  6.7).  It  is  emphasized,  however,  that  the 
algorithms  developed  under  this  contract  can  be  utilized  with  any  model 
capable  of  formulation  according  to  equations  (1)  through  (3). 


2.3.2  Polymeric  composites  at  low  homologous  temperatures  can  be  modeled  using 
internal  state  variable  theory  as  well.  However,  in  this  case  the  ISV's  are 
assumed  to  represent  locally  averaged  measures  of  various  damage  mechanisms 
such  as  matrix  cracking  and  interply  delamination.  The  constitutive  equations 
are  given  by  [ 15,16) 
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where  a^.  is  the  residual  stress  tensor,  is  the  elastic  modulus 

tensor,  I^i  the  damage  modulus  tensor  for  each  damage  mode,  and 
n  ranges  from  1  to  the  number  of  damage  modes.  For  example,  matrix  cracks, 
interply  delaminations,  and  fiber  breaks  each  represent  one  damage  mode. 
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The  internal  state  variables  are  described  by  history  dependent  ISV 
growth  laws  of  the  form 

*1j  =  Qij(£ka*  ekt*  T‘  “ka^ 

Equations  of  the  form  of  (4)  and  (5)  have  been  utilized  in  the  analysis 
of  composite  space  structures  in  Appendices  6.5,  6.8,  and  6.9. 

2.3.3  Metal  matrix  composites  are  expected  to  be  utilized  commonly  in  space 
structural  applications  due  to  their  high  melting  temperatures  compared  to 
polymeric  composites.  No  appropriate  constitutive  equations  for  these 
materials  were  found  in  the  literature.  It  was  therefore  felt  that  some 
constitutive  model  development  was  warranted  for  this  class  of  materials.  The 
distinguishing  feature  of  metal  matrix  composites  is  the  substantial  inelastic 
(either  elastic-plastic  or  viscoplastic)  nonlinearity  which  occurs  in  the 
matrix.  On  the  other  hand,  chopped  fiber  metal  matrix  composites  do  not 
exhibit  the  degree  of  layered  anisotropy  observed  in  laminated  continuous 
fiber  polymeric  composites.  Due  to  these  differences,  the  internal  state  in 
metal  matrix  composites  can  be  significantly  different  from  polymeric 
composites.  Accordingly,  a  generalized  model  was  developed  for  this 

material.  Although  the  model  is  an  extension  of  previous  research  on 
polymeric  composites  [151,  the  mechanics  of  damage  development  are  totally 
different.  The  details  of  this  model  are  given  in  references  17  and  18.  In 
addition,  a  synopsis  is  given  in  Appendix  6.6. 

The  constitutive  framework  is  based  on  a  continuum  mechanics  approach 
with  constraints  on  the  relations  provided  by  thermodynamics  and  fracture 
mechanics.  The  general  model  is  applicable  to  materials  with  damage  (such  as 


(10) 
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Therefore,  expanding  equations  (7)  and  (8)  in  Taylor  series  expansions  in 
terms  of  their  arguments,  substituting  into  equation  (10)  and  truncating 
higher  order  terms  results  in 


where  E|_  is  the  initial  loading  elastic  modulus.  Now  define  the  initial 
unloading  modulus  Ey  such  that  (See  Fig.  1): 
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It  is  assumed  that  at  relatively  low  homologous  temperatures  the 
inelastic  strain  remains  constant  when  unloading,  so  that: 
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Assuming  linear  elastic  unloading  of  the  matrix,  the  change  in  damage  is 
proportional  to  the  change  in  strain: 
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temperature  change. 
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—  =  constant  (upon  unloading)  =  6 

o  6 

Therefore,  for  unloading: 


Zu  =  El(1-b)  (16) 

The  model  described  above  provides  the  motivation  for  the  experimental 
research.  It  is  hypothesized  that  the  two  parameters  that  are  defined  in  this 
model  (a^.  and  e)  can  be  determined  by  experimental  methods.  By  determining 
the  change  between  the  initial  loading  and  unloading  moduli  of  the  composite 
in  a  uniaxial  mechanical  test,  s  (defined  in  equation  16)  can  be  found.  It 
is  also  observed  that  a.,  (defined  in  equation  9)  can  be  determined  by 

1  J 

evaluating  the  amount  of  surface  area  in  the  composite.  It  is  therefore 
desirable  to  determine  if  a  cause-and-effect  relationship  exists  between  the 
mlcrostructural  damage  (a..)  and  the  stiffness  loss  (b). 

*  J 

The  primary  objective  of  the  experimental  effort  was  to  develop  a 
technique  for  determining  and  evaluating  damage  in  metal  matrix  composites. 
This  technique  was  required  to  be  capable  of  detecting  cracks  and  voids  (free 
surfaces)  in  the  composite.  These  cracks  are  generally  on  the  order  of 
microns  in  characteristic  dimension,  so  that  scanning  electron  microscopy  was 
required  to  measure  the  damage.  Specimens  were  loaded  to  different  levels  and 
the  damage  was  studied  at  each  increment.  Once  the  amount  of  damage  was 
determined  it  could  be  input  into  the  general  constitutive  model  for  metal 
matrix  composites. 

The  material  used  in  this  study  was  obtained  from  ARCO  Metals  Silag 
Operation  in  Greer,  S.C.  The  composition  of  the  material  is  6061  Aluminum 
with  a  twenty  percent  volume  fraction  of  F-9  silicon  carbide  whiskers.  Plate 


was  made  from  the  materials  by  a  powder  metallurgy  process  and  cast  into 
billets.  The  billets  were  then  rolled,  extruded  or  machined  to  the  desired 
shapes.  The  SiC  whiskers  average  two  microns  in  diameter  and  twenty  microns 
long.  The  composite  has  a  T-6  temper.  Tensile  test  coupons  were  machined  in 
accordance  with  ASTM  E-8  (Tension  testing  of  metallic  materials)  to  the 
dimensions  shown  in  Fig.  2.  For  the  initial  portion  of  the  study  all 
specimens  were  machined  with  the  same  orientation  with  respect  to  the  plate 
for  the  purpose  of  uniformity  (with  respect  to  the  SiC  whisker  orientation). 
A  second  phase  of  the  testing  involve  the  use  of  tensile  test  specimens 
oriented  perpendicular  to  the  initial  specimens.  Details  of  the  experimental 
procedure  are  given  in  Appendix  6.6. 

It  was  shown  in  the  development  of  the  constitutive  model  that  the  damage 
parameter,  a.^,  can  be  used  to  predict  stiffness  losses.  Therefore,  it  was 
the  objective  or  the  experimental  research  to  measure  stiffness  loss  and  the 
associated  microstructural  damage  as  a  function  of  strain  level  in  order  to 
qualitatively  assess  the  applicability  of  the  model  to  the  Al-SiC  metal  matrix 
composite.  This  experimental  objective  was  carried  out  by  determining  the 
initial  loading  and  subsequent  unloading  moduli  of  tensile  specimens  oriented 
parallel  (0°)  and  perpendicular  (90°)  to  the  principle  rolling  direction  of  a 
plate  fabricated  from  6061-T6  aluminum  with  silicon  carbide  whiskers.  In 
addition,  scanning  electron  microscopy  was  utilized  to  characterize  and 
quantify  load-induced  changes  in  the  microstructural  damage  associated  with 
the  silicon  carbide  particles. 

The  results  showed  that  the  AT -Si C  plate  was  anisotropic  with 
approximately  15-20%  difference  in  the  moduli  of  the  specimens  oriented  in  the 
0°  and  90°  directions.  Also,  the  SEM  photomicrographs  indicated  that  the  SiC 
whiskers  were  oriented  more  or  less  parallel  to  the  principle  rolling 
direction. 
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There  was  very  little  difference  between  the  initial  loading  and 
subsequent  unloading  moduli  of  the  0°  specimens,  as  shown  in  Fig.  3.  Also, 
there  was  no  apparent  load-induced  change  in  the  state  of  microstructural 
damage  in  these  specimens.  On  the  other  hand,  there  was  a  significant 
reduction  in  the  moduli  of  specimens  oriented  in  the  90°  direction,  as  shown 
in  Fig.  4.  Furthermore,  the  photomicrographs  revealed  very  obvious  and 
significant  load-induced  changes  in  the  state  of  microstructural  damage  in  the 
90°  specimens. 

The  results  illustrate  a  clear  cause  and  effect  between  the  Increase  in 
load-induced  microstructural  damage  and  a  decrease  in  the  elastic  modulus  of 
the  Al-SiC  metal  matrix  composite.  It  is  concluded  from  these  results  that 
the  constitutive  behavior  of  a  short-fiber  reinforced  metal  matrix  composite 
can  only  be  modelled  by  an  appropriate  treatment  of  the  microstructural  damage 
associated  with  the  fiber  particles.  Although  the  model  developed  herein  is 
capable  of  accounting  for  these  effects,  due  to  the  qualitative  nature  of  the 
results  obtained,  the  constitutive  equations  for  metal  matrix  composites  have 
not  at  this  time  been  utilized  to  model  the  response  of  large  space 
structures. 

2.4  Coupled  Energy  Balance  Law 

The  energy  balance  law  for  thermomechanically  coupled  media  of  the  type 
described  in  Section  2.3.1  has  been  constructed  [19]  (see  also  Appendix 
6.1].  This  equation  can  be  utilized  to  predict  temperature  rise  in  a 
thermoviscoplastic  medium  subjected  to  cyclic  mechanical  loading.  This 
equation  Is  in  general  a  statement  of  conservation  of  energy  and  represents  a 
modification  of  the  Fourier  heat  conduction  equation  given  by 


Comparison  of  Initial  Loading  and  Subsequent  Unloading  Moduli  for  the 
0°  Specimen  Orientation 
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where  a-|kl  is  the  internal  state  variable  modeling  plastic  strain,  akji  is  the 
thermal  conductivity  tensor,  p  is  the  mass  density,  Cv  is  the  specific  heat, 
q^  is  the  heat  flux  vector,  and  r  is  the  specific  heat  supply. 

The  above  result  has  been  utilized  to  predict  the  thermomechanical 
response  of  a  single  perfectly  insulated  truss  element  to  cyclic  mechanical 
loading  (see  Appendix  6.1).  As  shown  in  Figs.  5  and  6,  substantial 
temperature  rise  (approximately  3.7°C)  is  predicted  for  each  cycle. 

In  polymeric  composites  the  majority  of  the  strain  energy  lost  to 
inelastic  deformations  may  be  expended  in  the  creation  of  internal  surfaces 
called  damage.  It  is  therefore  assumed  to  be  unnecessary  to  construct  two- 
way  coupled  energy  balance  laws  for  these  materials  and  the  classical  Fourier 
heat  conduction  equation  is  adequate  for  modeling  the  temperature  field. 
Therefore,  the  models  developed  herein  utilize  only  one  way  coupling  for 
polymeric  composite  media;  that  is,  the  temperature  field  affects  the 
displacement  field  but  not  vice  versa. 


2.5  Space  Structural  Response  Algorithms 

Due  to  the  nonlinearity  introduced  by  the  constitutive  equations 
developed  in  Section  2.3,  as  well  as  radiation  thermal  boundary  conditions, 
approximate  techniques  must  be  utilized  in  order  to  obtain  results  for 
geometries  representing  space  structures.  Accordingly,  the  following  finite 
element  computer  algorithms  were  developed  during  the  course  of  the  research 
effort: 
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1)  one-dimensional  code  for  analysis  of  two  way  coupled  thermoviscoplastic 
media  (see  Appendices  6.3  and  6.4); 

2)  truss  code  for  analysis  of  composite  space  structures  with  spatially 
variable  and  history  dependent  damage  (see  Appendix  6.5); 

3)  two-dimensional  continuum  code  for  predicting  quasi-static  response  of 
elastic-plastic  media  to  mechanical  inputs  (see  Appendix  6.2); 

4)  frame  code  with  one-way  coupled  thermal  analysis  for  predicting  response  of 
composite  space  structures  to  thermomechanical  inputs  (see  Appendix  6.8); 
and 

5)  beam  code  for  analysis  of  beam-like  composite  space  structures  to  spatially 
and  history  dependent  damage  (see  Appendix  6.9). 

Details  of  the  models  are  given  in  the  various  appendices  cited  above. 

In  addition,  further  information  on  the  algorithm  mentioned  in  item  4)  above 

is  contained  in  reference  20. 


2.6  Model  Results  for  Large  Soace  Structures 


For  purposes  of  illustrating  the  capability  of  the  models  developed  under 
the  contract,  several  sample  problems  are  provided  here  for  representative 
space  structures.  These  examples  fall  into  the  following  four  categories:  I) 
heat  generation  due  to  cyclic  loading  of  metallic  members;  2)  frequency  and 
mode  shape  degradation  of  composite  truss  structures;  3)  radiation  induced 
response  of  composite  frame  structures;  and  4)  degradation  of  dynamic  response 
of  composite  beams  with  damage.  Results  for  these  four  examples  are  discussed 
briefly  below. 


2.6.1  Heat  generation  in  metallic  members  occurs  due  to  coupling  between 
thermal  and  mechanical  effects,  as  discussed  in  Section  2.4.  Considerable 
research  was  performed  on  this  subject,  as  detailed  in  Appendices  6.1,  6.3, 
6.4,  and  6.7.  These  results  are  summarized  here. 


It  was  found  that  for  an  insulated  truss  member  composed  of  a 
representative  metal,  when  the  member  is  intentionally  loaded  cyclically  to 
the  postyielded  state  (in  order  to  induce  significant  material  damping)  shown 
in  Fig.  5,  a  temperature  rise  of  3.7°K  occurred  on  each  load  cycle,  as  shown 
in  Fig.  6.  This  was  cause  for  further  study,  since  a  modal  response  of  the 
structure  could  eventually  result  in  structural  component  melt-down. 
Therefore,  a  subsequent  effort  was  made  to  account  for  more  realistic  thermal 
boundary  conditions.  This  required  that  spatial  variability  of  the  field 
parameters  be  incorporated  into  the  model,  so  that  it  became  necessary  to 
utilize  the  finite  element  method.  This  resulted  in  a  highly  complex 
nonlinear  and  numerically  stiff  algorithm  due  to  the  viscoplastic  constitutive 
equations,  as  described  in  Appendix  6.3.  Results  indicated  that  for  the  case 
of  a  truss  member  with  insulated  longitudinal  boundaries  and  nonzero  thermal 
flux  at  each  end,  for  the  loading  input  shown  in  Fig.  5  the  temperature  rise 
per  cycle  was  reduced  to  1.0°K  per  cycle,  as  shown  in  Fig.  7.  Since  this  was 
still  considered  to  be  sufficient  to  lead  to  structural  failure,  it  was 
decided  to  incorporate  the  effect  of  nonlinear  radiation  boundary  conditions 
to  the  algorithm.  As  shown  in  Figs.  8  and  9,  the  temperature  rise  is  quite 
substantial  for  fifty  cycles,  even  at  moderate  frequencies  (cases  I  and  II 
represent  different  member  coatings),  reaching  90°K  for  a  frequency  of  5  Hz. 
Finally,  bounds  on  the  predicted  temperature  rise  versus  stress  amplitude  are 
shown  in  Fig.  10.  Since  this  amount  of  heating  is  unacceptable  in  metallic 
structures,  it  is  concluded  that  the  intentional  use  of  metal  inelasticity  to 
induce  passive  structural  damping  can  possibly  lead  to  catastrophic  heating  of 
the  structure.  Therefore,  attempts  to  induce  damping  via  this  mechanism 
should  be  viewed  cautiously. 
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2.6.2  Degradation  of  modal  frequencies  and  shapes  can  occur  in  composite  space 
structures  wherein  microstructural  damage  occurs  as  a  result  of  the  history  of 
loading.  As  an  example,  consider  the  truss  structure  shown  in  Fig.  11.  This 
beam-like  structure  is  cantilevered  at  one  end  to  simulate  an  antenna  boom. 
The  structure  is  sixty  feet  long  with  bays  ten  feet  long  by  three  feet  wide. 
The  structure  is  constructed  from  graphite-epoxy  composite  material  with  a 
quasi-isotropic  ply  layup.  Experimental  research  indicates  that  the  material 
may  undergo  up  to  15  percent  loss  in  stiffness  due  to  cyclic  thermomechanical 
fatigue  which  causes  a  variety  of  damage  modes  in  the  structure.  Additional 
loss  of  stiffness  may  be  attributed  to  elevated  temperature  and  chemical 
changes  due  to  solar  radiation  and  other  environmental  efffects.  In  this 
model  the  properties  are  degraded  spatially  on  an  element  by  element  basis  as 
a  function  of  the  stress  history  in  the  structure  induced  by  long  term 
thermomechanical  cyclic  loading.  Stress  amplitudes  were  obtained  by  using 
displacements  corresponding  to  the  first  modal  shape  and  the  degraded 
properties  were  computed  by  assuming  a  linear  damage  law  bases  on  peak  stress 
amplitude.  Because  the  boom  is  fixed  on  one  end,  the  stresses  are  highest 
there  and  stiffness  degrades  the  most  at  the  fixed  end.  Modal  frequencies  and 
shapes  were  then  computed  for  the  five  cases  where  the  maximum  degradation 
within  the  structure  was  5%,  10%,  ...,  and  25%.  Fig.  12  indicates  the 
decrease  in  the  fundamental  frequencies  of  the  first  two  resonant  modes  as  a 
function  of  this  spatially  induced  damage.  Fig.  13  shows  that  the  shape  of 
the  first  mode  undergoes  no  appreciable  change  in  shape  as  the  damage  occurs, 
which  is  due  to  the  fact  that  the  first  mode  is  a  symmetric  mode.  The  second 
mode  shape,  however,  shown  in  Fig.  14,  undergoes  a  substantial  change  in 
shape.  This  result  indicates  that  active  control  mechanisms  which  may  be 
placed  according  to  the  original  undamaged  mode  shapes  may  not  be  capable  of 
controlling  all  modes  as  the  dynamic  response  of  the  structure  changes  over 
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Fig.  13  Change  in  First  Mode  Shape  of  Typical  Space  Truss  Structure  Due  to 
Material  Property  Degradation. 
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Fig.  14  Change  in  Second  Mode  Shape  of  Typical  Space  Truss  Structure  Due  to 
Material  Property  Degradation. 


long  time  periods.  The  control  designer  must  be  cognizant  of  these  modal 
changes  if  he  is  to  design  a  workable  control  system.  Adaptive  or  robust 
control  systems  will  be  required. 

Numerous  other  results  are  presented  in  Appendix  6.5.  For  example,  it  is 
shown  that  if  the  material  properties  degrade  on  only  one  side  of  the 
structure,  as  might  happen  due  to  solar  radiation,  the  mode  shapes  are 
completely  changed. 

In  conclusion,  these  preliminary  results  seem  to  indicate  that  small 
changes  (or  errors)  in  material  properties  as  they  change  or  degrade  due  to 
fatigue  damage,  etc.  may  produce  significant  changes  in  predicted  frequency 
and  modal  response.  Correspondingly,  this  affects  the  ability  to  design 
effective  control  systems  and  places  an  even  greater  burden  on  the  control 
designer  to  develop  systems  which  account  for  these  structural  changes.  It  is 
clear  that  an  understanding  of  material  behavior  in  space  environments  and  its 
impact  on  structural  response  is  very  important  to  successful  design  and 
development  of  large  space  structures. 

2.6.3  Radiation  Induced  response  of  composite  frame  structures  is  caused  by 
thermal  strains  resulting  from  solar  and  earth  radiation.  The  thermoelastic 
boundary  value  problem  is  complicated  by  several  factors.  First,  a  one  way 
coupling  between  temperature  and  displacements  exists.  It  is  one  way  coupled 
in  that  displacements  depend  on  temperatures.  Secondly,  the  problem  is 
nonlinear  due  to  the  introduction  of  radiation  boundary  conditions.  Thirdly, 
there  are  constantly  changing  thermal  loading  conditions  due  to  varying  earth- 
structure-sun  orientation.  Finally,  geometrical  factors  such  as  shadowing  and 
interelement  radiation  and  conduction  exist.  These  factors  combine  to  create 
a  highly  complex  problem,  as  shown  in  Fig.  15. 
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Fig.  15  A  Structural  Element  in  Space  Environment 


Past  research  indicates  that  a  structural  member,  modeled  as  a  thin- 
walled  graphite/epoxy  tube  under  a  constant  solar  flux,  experiences  a 
significant  temperature  gradient  through  its  cross-section  due  to  the  low 
thermal  conductivity  of  the  material.  This  gradient  leads  to  bending  in  the 
structural  member  itself,  a  condition  normally  neglected  in  thermal/structural 
analysis  models.  This  response  is  important  for  two  reasons.  First,  the 
bending  of  a  structural  member  reduces  the  maximum  buckling  load  that  member 
can  sustain.  Second,  thermally  induced  vibrations  may  lead  to  fatigue,  which 
is  important  in  predicting  the  long-term  behavior  of  the  material. 

The  purpose  of  this  research  was  the  development  of  an  integrated,  one¬ 
way  coupled  thermoelastic  model  for  transient  analysis  of  large  composite 
space  structures.  The  primary  load  source  is  thermal  strains  induced  by  solar 
and  earth  radiation.  Therefore,  the  model  must  be  capable  of  transforming 
these  thermal  loads  into  their  mechanical  equivalents.  Due  to  the  presence  of 
radiation  and  temperature  dependent  material  properties,  the  model  is  highly 
nonl inear. 

The  thermoelastic  analysis  is  a  completely  numerical  one.  The  model 
satisfies  all  the  requirements  for  a  fully  integrated  thermal/structural 
analysis  model.  These  are:  a  common  finite  element  methodology,  utilization  of 
a  single  geometric  model,  improved  thermal  analysis,  minimized  data  transfer 
between  analyses,  and  a  thermal  analysis  fully  integrated  into  the  structural 
analysis.  However,  the  model  is  unique  in  two  key  areas.  First,  temperature 
gradients  across  member  cross-sections  are  accounted  for.  Thermal  moments  and 
extensions  are  calculated  directly  form  integration  of  the  resulting 
temperature  fields.  Second,  structural  members  are  modeled  with  beam  elements 
enabling  bending  in  the  structural  members  themselves. 
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In  performing  the  thermoelastic  analysis,  it  is  assumed  that  temperatures 
and  displacements  are  one-way  coupled.  That  is,  temperature  fields  within 
structural  members  may  be  determined  independently,  then  used  as  input  to  the 
structural  analysis.  The  algorithm  is  as  follows.  For  a  given  time  step,  the 
proper  heat  loads  are  evaluated.  Then,  finite  elements  are  used  to  construct 
the  temperature  field  through  member  cross-sections.  A  standard  two- 
dimensional  finite  element  formulation  is  applied  to  the  following  equations 
governing  heat  transfer  in  the  domain  of  a  cross-section,  (see  Fig.  16): 
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Here  A  is  cross-sectional  area,  B  is  the  boundary,  t  is  time,  p  is  material 
density,  Cv  is  the  specific  heat,  kx  and  ky  are  the  thermal  conductivities,  T 
is  temperature,  Q  is  internal  heat  generation,  nx  and  ny  are  vector  normals,  q 
is  the  flux,  h  is  the  film  coefficient,  e  is  emissivity,  S  is  Bolt2man's 
constant,  and  Ta  and  Tr  are  ambient  and  reference  temperatures.  The  resulting 
temperature  field  is  then  converted  into  bending  moments  and  extensions  using 
the  following  equations: 
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Fig.  16  Heat  Transfer  in  a  Selected  Cross-Section 


the  coefficient  of  thermal  expansion,  T  is  the  current  temperature  field,  and 
Tq  is  the  initial  temperature  field  or  the  temperature  field  in  the 
unstrained  state.  The  thermal  force  and  moments  about  the  local  y  and  z  axes 
are  given  by  P-j-,  M  and  ,  respectively.  These  loads  are  initially 
calculated  in  the  local  coordinates  of  their  cross-section,  then  transformed 
into  global  coordinates  to  be  used  as  forcing  functions  in  the  structural 
analysis.  The  structural  analysis  portion  of  the  model  is  for  linear  space 
frame  geometries.  It  is  the  result  of  applying  the  standard  finite  element 
formulation  to  the  governing  differential  equations  of  beam  motion: 


Mu,t),t  -  (EAU,X),X  =  PT,XX 

(23) 
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(26) 

Here  u,v  and  w  are  displacements  in  the  x,y,  and  z  coordinate  directions,  e 
is  the  rotation  about  the  x  axis,  I  is  the  mass  inertia  of  a  cross-section, 
G  is  the  modulus  of  rigidity,  J  is  the  polar  moment  of  inertia,  1^  and 
are  the  bending  inertias  about  the  y  and  z  axes,  and  PT,  and  M^T,  and  Mz  are 
given  in  equations  (20),  (21),  and  (22).  Once  the  deformations  and  stresses 
are  determined,  time  is  incremented  and  the  process  repeated.  Forward 
integration  in  time,  within  the  structural  analysis,  is  via  the  Newmark 
method.  The  thermal  analysis  utilizes  the  Crank-Nichol son  method.  The  final 
algorithm  is  described  schematically  in  Fig.  17. 
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The  boom  shown  in  Fig.  11  was  analyzed  for  two  physical  cases.  In  case 
one,  the  structure  was  assumed  to  be  in  thermal  equilibrium  and  stress-free  in 
sunlight.  At  time  t=0,  the  structures  moves  into  shadow.  In  case  two,  the 
opposite  occurrs,  that  is,  the  unstressed  structure  moves  from  shadow  to 
sunlight.  The  maximum  axial  stress  divided  by  the  yield  stress  is  plotted 
against  time  for  a  typical  member  in  Figs.  18  and  19.  It  can  be  seen  from  the 
results  that  bending  stresses  are  not  negligible,  indicating  that  the 
structure'  must  be  treated  as  beam-like  rather  than  truss-like  in  nature. 
Furthermore,  the  long  time  necessary  to  induce  significant  stresses  indicates 
that  even  though  the  transition  from  shadow  to  sunlight  is  assumed  to  be 
instantaneous,  no  inertial  effects  are  induced  by  radiation.  Further  results 
from  this  portion  of  the  research  are  given  in  Appendix  6.8. 


2.6.4  Damage  in  composite  beams  can  substantially  alter  their  dynamic 
response.  Fibrous  composites  are  known  to  undergo  a  small  but  significant 
amount  of  stiffness  loss  due  to  load  induced  microcracking.  This  stiffness 
loss  usually  occurs  over  several  hundred  thousand  load  cycles.  Due  to  the 
stress  dependent  nature  of  the  damage,  the  stiffness  loss  is  spatially 
variable  and  concentrated  in  the  areas  of  high  stresses.  This  spatial  change 
in  the  material  properties  of  the  structure  results  in  appreciable  changes  in 
the  dynamic  response  of  the  structure. 

Although  the  authors  have  previously  developed  a  qualitative  method  for 
predicting  response  of  structures  of  this  type  (see  appendix  6.5),  it  is  not 
possible  to  construct  a  precise  history  dependent  structural  algorithm  for 
complex  space  structures  due  to  the  excessive  computational  times  required  to 
obtain  accurate  results.  A  more  accurate  method  has  now  been  developed  for  a 
single  beam  member  with  various  boundary  conditions.  Since  this  model  carries 


structures,  it  is  possible  to  determine  the  actual  structural  response  for  a 
load  input  of  several  hundred  thousand  cycles.  the  following  is  a  brief 
description  of  this  procedure. 

The  well  known  partial  differential  equation  for  the  free  vibration  of  a 
beam  is 

— 2  (El  £$)  +  oA  =  0  (27) 

3X  3X  3t 

where  E  is  Young's  modulus,  I  is  the  moment  of  inertia  of  the  cross-section,  A 
is  the  cross-sectional  area,  p  is  the  mass  density,  y  is  the  transverse 
displacement,  x  is  the  axial  coordinate,  and  t  is  time. 

A  number  of  solutions  to  the  above  differential  equation  are  available  in 
the  literature  for  both  uniform  (constant  cross-section)  and  nonuniform 
(variable  cross-section)  with  different  boundary  conditions.  Most  of  the 
solutions  are  for  beams  with  homogeneous  material  properties.  These  solutions 
have  been  obtained  by  assuming  that  the  stiffness  of  a  structural  element  is 
constant  in  time  and  therefore  independent  of  loading  history.  Neither 
material  damage  nor  environmental ly  caused  degradation  are  considered  in  these 
analyses. 

Due  to  the  occurrence  of  load  induced  and  history  dependent  damage  in 
composite  mateials,  these  previously  obtained  results  represent  unrealistic 
approximations  of  the  actual  structural  behavior.  In  particular,  the  resonant 
frequencies  and  mode  shapes  of  the  structure  can  be  severely  altered  by  the 
introduction  of  spatially  variable  damage  and  environmentally  caused 
degradation,  the  stiffness  of  a  structure  is  no  longer  a  constant,  since  it 
will  change  substantially  according  to  the  stress  distribution  and  the  history 
of  external  loading.  The  stiffness  loss  may  change  the  natural  frequencies 


and  mode  shapes  substantially.  With  the  material  damage  and  environmental ly 
caused  degradation  involved,  the  differential  equation  becomes  difficult  if 
not  impossible  to  solve  in  closed  form. 

The  concept  of  Internal  state  variables  (ISV)  is  introduced  to  represent 
the  history  dependent  change  of  stiffness.  An  internal  state  variable  D  is 
utilized  as  a  local  ISV  representing  the  damage  state  together  with  the  ISV 
growth  law,  the  finite  element  solution  technique  can  be  modified  to  account 
for  the  history  dependent  stiffness  of  the  beam  element,  with  resulting  field 
equations 

(Ml  {y}  +  [ K ] {y }  =  {0}  (28) 

where 

M.  j  =qJL  N^jpAMdx 

K.j  =  qJL  E(e,T,0)I(x)N»i,jdx  (29) 

The  above  set  of  second  order  ordinary  differential  equations  for  each  element 
is  combined  to  represent  the  eigenvalue  problem  for  the  beam  structure. 

The  occurrence  of  damage  will  cause  the  loss  of  stiffness,  that  is  the 
stiffness  is  history  dependent.  Experimental  results  indicate  that  the  time 
scale  of  damage  and  degradation  is  very  long  compared  to  the  first  fundamental 
frequency  of  the  structure.  Therefore,  the  mathematical  algorithm  is  treated 
as  linear  with  slowly  varying  coefficients.  In  this  research,  particular 
interest  was  placed  on  the  natural  vibration  solution  of  a  beam  structure 
with  history  dependent  stiffness  and  the  investigation  of  the  possible  effect 
of  material  damage  and  stiffness  reduction  on  the  natural  frequencies  and  mode 
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shapes  of  planar  beam  structures  with  various  boundary  conditions  (free-free, 
clamped-free,  clamped-clamped  and  simply  supported). 

The  research  also  focused  on  the  investigation  of  the  internal  state 
variable  representation  of  the  damage  phenomenon.  The  damage  in  a  composite 
material  includes  a  sequence  of  microstructural  and  macrostructural  events 
such  a  microvid  growth,  matrix  cracking,  edge  delamination  and  fiber 
fracture.  The  most  significant  effect  of  damage  on  the  material  properties  is 
that  the  stiffness  will  be  substantially  changed  during  the  life  of  the 
component.  The  constitutive  equation  for  a  composite  material  could  be 
represented  as 


o  =  E(e  -  e  ) 


where  E  is  Young's  modulus,  which  will  change  according  to  the  damage  D  as 


E  =  E0(l  -  D) 


where  the  subscript  o  represents  the  initial  condition.  The  damage  D  is  an 
internal  state  variable  describing  the  damage  phenomenon  during  the  life  of 
the  composite  structure,  which  is  governed  by  the  internal  state  variable 
growth  law 


D  =  f(e,  T,  D) 


A  typical  result  for  the  degradation  if  the  frequency  of  the  first  mode 
of  the  simply  supported  composite  beam  shown  in  Fig.  20  is  shown  in  Fig.  21. 
Further  results  are  shown  in  Appendix  6.9  However,  it  was  found  that  the 
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mode  shapes  are  not  substantially  affected  by  damage.  We  believe  that  this 
was  due  to  the  assumption  that  the  damage  was  assumed  to  vary  through  the 
thickness  due  to  bending  stresses.  An  experimental  effort  is  currently 
underway  to  determine  the  accuracy  of  this  assumption. 

2.7  Conclusions 

This  research  effort  has  resulted  in  a  number  of  important  conclusions. 
We  summarize  a  few  of  these  as  follows: 

a)  in  metallic  structures  wherein  material  inelasticity  is  utilized  to  produce 
passive  damping,  thermomechanical  coupling  can  lead  to  castastrophic 
heating  of  the  structure; 

b)  moderate  changes  in  stiffness  of  composite  structural  members  can  lead  to 
substantial  changes  in  frequencies  and  mode  shapes  of  the  structure; 

c)  solar  and  earth  radiation  in  composite  structures  can  lead  to  thermal 
gradients  which  lead  to  substantial  bending  in  large  space  structures,  thus 
negating  the  efficacy  of  truss  analysis;  and 

d)  in  metal  matrix  composite  members  there  is  a  qualitative  comparison  between 
microstructural  damage  and  macrostructural  stiffness  loss. 

It  is  the  general  conclusion  of  these  researchers  that  material 
inelasticity  and  temperature  effects  cannot  be  disregarded  in  modelling  the 
dynamic  response  of  large  space  structures. 
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1986. 

5.2  Awards  and  Achievements 

1.  Dr.  Allen  has  been  named  Associate  Editor  of  the  Journal  of  Spacecraft 
and  Rockets. 

2.  Dr.  W.E.  Haisler  has  been  named  Head  of  the  Aerospace  Engineering 
Department  at  Texas  A&M  University. 

3.  The  textbook  entitled  Introduction  to  Aerospace  Structural  Analysis, 
co-authored  by  Drs.  Allen  and  Haisler,  has  been  published  by  John 
Wiley. 

4.  Dr.  W.E.  Haisler  has  been  named  to  the  Halliburton  Chair  at  Texas  A&M 
University. 

5.  Drs.  Allen  and  Haisler  have  been  named  Texas  Engineering  Experiment 
Station  Research  Fellows  for  1984-1985. 

6.  Dr.  Allen  has  received  the  General  Dynamics  Award  for  Outstanding 
Teaching  and  Research  in  the  College  of  Engineering  at  Texas  A&M 
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Professor. 

8.  Dr.  Allen  has  been  named  Texas  Engineering  Experiment  Station  Research 
Fellow  for  1985-1986. 
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Abstract 

A  thermodynamic  model  is  presented  for  predicting  the  thenaomechanical 
response,  including  temperature  change,  in  a  uniaxial  bar  composed  of  a 
thermoviscoplastic  metallic  medium.  The  model  is  constructed  using  the 
concept  of  internal  state  variables,  and  it  is  shown  that  this  general  frame¬ 
work  is  capable  of  encompassing  several  constitutive  models  currently  used 
to  predict  the  response  of  rate  sensitive  metals  in  the  inelastic  range. 

Results  are  obtained  for  monotonic  loading  which  agree  with  predicted  re¬ 
sults  previously  obtained  by  Cernocky  and  Krempl  for  mild  steel  at  room 
temperature.  The  model  is  then  utilized  in  conjunction  with  Bodner  and 
Partora's  constitutive  equations  to  predict  temperature  change  in  Inconel 
(IN)  100  subjected  to  both  monotonic  and  cyclic  loading  at  1005oK(1350oF) . 

Introduction 

It  has  long  been  known  that  mechanical  and  thermodynamic  coupling  ex¬ 
ists  in  solid  bodies  [1,2].  However,  in  elastic  bodies  this  coupling  is 
negligible  except  when  mass  inertia  is  not  negligible  due  to  flux  of  heat  generated 
through  the  boundary  of  the  body  [3].  However,  in  thermoviscoplastic  metals 
the  conversion  of  mechanical  energy  to  heat  may  be  singificant  even  under 
non-inertial  conditions,  especially  since  material  properties  become  ex¬ 
tremely  temperature  sensitive  in  the  inelastic  range  of  response  [4-11]. 

Similar  research  has  been  performed  on  non-metallic  media  [12-15]. 


General  continuum  mechanics  models  have  been  formulated  for  broad  classes 


of  materials  [16-19].  However,  to  this  author's  knowledge  only  recently 
has  attention  been  paid  to  the  coupled  heat  conduction  equation  for  thermo- 
viscoplastic  metals  [11,20].  Recently,  Cernocky  and  Krempl  [11]  have  pro¬ 
posed  a  model  for  predicting  the  temperature  rise  in  a  class  of  thermovis¬ 
coplastic  metals,  with  special  emphasis  on  test  coupons  subjected  to  either 
homogeneous  uniaxial  or  torsion  loadings.  In  this  paper  an  alternative  ap¬ 
proach  to  that  proposed  in  [11]  is  discussed.  This  method  uses  the  thermo¬ 
dynamics  with  internal  state  variables  originally  reported  in  [17]  and  dis¬ 
cussed  elsewhere  in  detail  for  metals  [18,21,22],  with  development  of  the 
multidimensional  coupled  heat-conduction  equation  in  [20]. 

The  research  herein  is  presented  in  three  parts:  field  formulation 
in  one-dimensional  form;  development  of  the  governing  equations  from  addi¬ 
tional  constitutive  assumptions;  and  numerical  results  for  selected  problems. 

Thermodynamics  of  a  Uniaxial  Thermoviscoplastic  Bar 

Consider  a  slender  bar  which  is  subjected  to  a  homogeneously  applied 
deformation  field  such  that  the  resulting  stress  field  is  everywhere  uniax¬ 
ial  in  the  x^=x  coordinate  direction,  as  shown  in  Fig.  1.  Rigor  would  re¬ 
quire  that  the  possibility  of  finite  deformations  be  considered.  However 
this  condition  is  covered  in  detail  elsewhere  [17,18,20,21,22],  and  for 
purposes  of  simplicity  only  infinitesimal  deformations  will  be  considered 
herein.  For  notational  simplicity,  then,  the  observable  mechanical  state 
variables  are 

u  =  u^  «  deformation  field,  (1) 

£  =  •  infinitesimal  strain  field,  and  (2) 


7  =  0,, 


stress  field 


(3) 


Although  transverse  components  of  deformation  and  strain  may  occur,  it  is 
assumed  that  they  are  not  necessary  tc  characterize  the  uniaxial  stress  a. 
The  mechanical  state  variables  (1)  through  (3)  are  adjoined  with  the 


thermodynamic  state  variables: 

e  E  internal  energy  per  unit  mass;  (4) 
r  E  heat  supply  per  unit  mass;  (5) 
s  E  entropy  per  unit  mass;  (6) 
T  E  absolute  temperature;  and  (7) 
q  E  q^  *  heat  flux  in  the  coordinate  direction  ,  (8) 


where  it  is  assumed  in  (8)  that  the  bar  is  isotropic  and  long  and  slender 
with  perfectly  longitudinal  boundaries  so  that  the  heat  flux  is  one-dimen¬ 
sional. 

In  accordance  with  the  theory  of  internal  state  variables  [17],  ob¬ 
servable  state  variables  (1)  through  (8)  are  now  supplemented  with  inter¬ 
nal  state  variable  growth  laws  in  order  to  characterize  the  state  of  inel¬ 
astic  bodies: 

E  scalar  valued  internal  state  variables,  k  *  1,  ...  ,  n;  (9) 
where  n  is  the  number  of  internal  state  variables  required  to  characterize 
the  state  of  the  body.  The  precise  nature  of  (9)  will  be  discussed  later. 

Parameters  (1)  through  (9)  are  assumed  to  be  functions  of  space  (x) 


c)  the  balance  of  energy, 

oe  -  oe  +  •  pr 

ox 

where  P  represents  the  mass  density;  and 

d)  the  second  law  of  thermodynamics. 


(12) 


(13) 


where  y  is  called  the  internal  entropy  rate  per  unit  mass. 

As  detailed  by  Coleman  &  Gurtin  [17],  equations  (10)  through  (13)  are 
now  supplemented  with  the  following  constitutive  assumptions: 


a  -  o(e,  T,  3X/3x,  a^)  ;  (14) 

e  ■  e(£,  X,  3T/3x,  a^)  ;  (15) 

s  -  s(e,  X,  3X73x,  a^)  ;  (16) 

q  «  q(e,  X,  3X73x,  a^)  ;  and  (17) 

-  a^e,  T*  3X/3x,  a^)  ,  (18) 


where  it  is  obvious  that  equations  (14)  through  (18)  satisfy  the  principle 
of  equipresence  [23].  Equations  (10)  through  (12)  and  (14)  through  (18) 
describe  eight  +n  equations  in  the  eight  +n  field  variables  u,  z,  a,  e, 
r,  s,  X,  q,  and  described  in  (1)  through  (9).  Xhese  are  adjoined  with 
boundary  conditions  on  the  surfaces  x  ■  0  and  x  *  L  to  prescribe  the  one¬ 
dimensional  field  problem. 

As  detailed  elsewhere  [17,18,20-22],  the  second  law  of  thermodynamics 
[inequality  (13)]  will  constrain  constitutive  assumptions  (14)  through  (18). 
Xhis  is  accomplished  by  defining  the  Helmholtz  free  energy: 

h  =  h(e,  X,  3 Xy3x,  c^)  -  e  -  Xs  ->  e  -  h  +  Is  ,  (19) 


which  together  with  the  Clausuis-Duhem  inequality  will  lead  to  the  conclu¬ 


sions  that 


h  -  h(e,  T,  ct^) 


Tf  *s<e-  T-  \> 


a  *  ^  3e  ,  and 


q  »  -k 


where  k  is  the  coefficient  of  heat  conduction  in  the  x^  coordinate  direc¬ 


tion.  Therefore,  equations  (19)  through  (23)  replace  equations  (14)  through 


(18)  as  more  concise  statements  of  the  constitutive  behavior,  and  it  can 


be  seen  that  specification  of  the  Helmholtz  free  energy  will  complete  the 


description  of  the  field  problem. 


Combination  of  equations  (12)  and  (19)  through  (23)  will  result  in  the 


coupled  heat  conduction  equation: 


3h  •  _  3*h  •  _  3*h  •  _  32 

0  3cl  “k  "  pT  3a,  3T  “k  "  pT  3e3T  e  ”  DT  3T 


3T7  1  +  3x 


where  summation  on  the  range  of  k  is  implied. 


Henceforth  in  this  investigation  it  will  be  assumed  that  there  is 


no  internal  heat  source  (other  than  material  dissipation)  so  that  r  =  0 


in  equation  (24).  In  addition,  it  will  be  assumed  that 


boundary  conditions  are  applied  in  such  a  way 


that  heat  flux  is  negligibly  small  and  q  *  0  in  equation  (24) .  This  last 


assumption  is  not  valid  under  most  physical  circumstances.  However,  it 


can  be  said  that  on  the  basis  of  heat  conduction  equation  (24)  neglecting 


heat  flux  will  result  in  an  upper  bound  for  the  temperature  rise  during 


mechanically  induced  energy  dissipation.  Inclusion  of  this  term  results 


in  a  spacially  dependent  boundary  value  problem  which  is  beyond  the  scope 


of  the  current  research.  However,  the  one-dimensional  model  proposed  herein 


does  encompass  the  heat  flux  phenomenon,  and,  as  such,  will  be  the  subject 
of  a  future  paper  by  the  author. 

Development  of  Governing  Equations  from 
Additional  Constitutive  Assumptions 


In  order  to  construct  the  Helmholtz  free  energy  function  the  elastic 
strain  is  first  defined  to  be 

E  — 

t  =  e  -  -  otS  ,  (25) 

where  is  the  total  inelastic  strain  in  the  x^  coordinate  system  [24], 

a  is  the  coefficient  of  thermal  expansion  in  the  x^  coordinate  direction, 

and  0  =  T  -  T  ,  where  T  is  the  initial  temperature  at  which  no  strain  is  ob- 
K  R 

served  under  zero  mechanical  load.  The  inelastic  strain  will  be  discussed 
in  greater  detail  in  the  next  section. 

It  is  now  postulated  that  the  Helmholtz  free  energy  may  be  expanded 
about  the  initial  configuration  in  terms  of  elastic  strain  and  temperature 
as  follows: 


h  -  >,  +  pe2  Cv  q2 

h  "  \  +  TB  E  'It  9 

where  the  subscript  R  denotes  the  equilibrium  value,  and 
h^  =  free  energy  in  state  R  «  constant. 


(26) 

(27) 


E  =  Young's  modulus  in  the  x^  coordinate  system,  (28) 

Cy  =  -TO^/ST2)  ■  specific  heat  at  constant  volume.  (29) 

E 

Mote  that  although  the  first  order  terms  in  e  and  6  have  been  neglected 


due  to  symmetry  conditions  due  to  the  form  of  equation  (25)  coupling  is 


retained  between  total  strain,  inelastic  strain  and  temperature.  Note  also 


that  the  energy  dissipation  due  to  microstructural  change  has  been  neglected 
in  free  energy  equation  (26)  because  this  mechanism  has  been  shown  to  contribute 


•  .N- v -.v 


7 


only  a  small  portion  of  energy  (<10%)  to  the  dissipation  process  [25]. 

Further,  the  fracture  energy  loss  due  to  microvoid  growth,  grain  boundary 
sliding,  and  intergranular  macrofracture  is  neglected  due  to  the  small 
strains  considered  herein. 

Although  the  second  order  Taylor  series  expansion  of  the  Helmholtz 
free  energy  given  in  equation  (26)  may  not  be  adequate  for  characterizing 
the  response  of  many  materials,  it  will  be  shown  in  the  next  section  that 
the  above  equations  are  a  suitable  framework  for  describing  the  material 
behavior  of  the  class  of  materials  considered  herein. 

Substitution  of  equation  (26)  into  energy  balance  law  (24)  and  uti¬ 
lizing  equation  (25)  will  result  in  the  coupled  heat  equation: 

[(Ee  -  Eai  +  Eal^o^  +  Ea2TT]  -  EaTe.  -  oCvT  «  0  ,  (30) 

where  the  terms  in  brackets  arise  due  to  inelastic  response  and  the  following 
term  is  the  classical  elastic  coupling  term  [3].  Equation  (30)  may  be  written 
in  the  following  equivalent  form: 


(Ee  -  Ea^  +  EaT^)a^  -  EaTe 

PC  -  Ea2T 
v 


(31) 


In  order  to  obtain  the  stress-strain  relation  the  Helmholtz  free  energy 
equation  (26)  may  be  substituted  into  equation  (22)  to  obtain 

a  -  E(e  -  ctj  -  a8)  .  (32) 


Equations  (31)  and  (32)  together  with  internal  state  variable  growth 
laws  (18)  will  be  sufficient  to  characterize  the  response  of  the  uniaxial 
bar  subjected  to  uniaxial  homogeneous  mechanical  loading  considered  herein. 


Selected  Problems  and  Numerical  Results 


It  has  been  shown  that  stress-strain  relation  (32)  together  with  in¬ 


ternal  state  variable  growth  laws  (18)  are  equivalent  to  several  models 


J .,  MM  ■>  ■>  |T«,  ■  ,  .,., ,v  ,  t* T. vny.  V.’1* v"HP»yw\rTVT 


recently  proposed  for  thermoviscoplastic  metals  [24].  These  include  Cer¬ 
nocky  and  Krempl  [11,26],  Valanis  [27],  Krieg,  et  al.  [28],  and  Allen  and 
Haisler  [29].  It  can  also  be  shown  that  several  others  are  in  accordance 
with  the  model  developed  herein  [30-34].  To  illustrate  this  point  two 
models  have  been  selected  for  further  discussion. 

Cernocky  and  Krempl* s  stress-strain  relation  may  be  written  in  the 
following  uniaxial  form: 

a  +  K(a,  e,  T)a  -  G(e,  t)  +  M(a,  e,  t)[c  -  ai]  ,  (33) 

where 

E  =  M/K  ,  (34) 

and  parentheses  imply  dependence  on  the  current  values  of  the  quantities 
enclosed.  Equations  (33)  and  (34)  can  be  shown  to  be  in  agreement  with 
stress-strain  equation  (32)  by  defining  the  inelastic  strain  such  that 

ax  -  [a  -  G(e, T )]/M(a,  e, T )  ,  (35) 

so  that 

t 

oi^t)  -  /  OjCt^dt’  ,  (36) 

CR 

where  t^  is  the  reference  time,  t'  is  a  dummy  variable  of  integration,  and 
t  is  the  time  of  interest.  Thus,  since  G,  K  and  M  are  not  history  depen¬ 
dent,  Cernocky  and  Krempl's  model  is  a  single  internal  state  variable  model 
and  equations  (31),  (32),  (35)  and  (36)  describe  the  uniaxial  bar  problem 
using  Cernocky  and  Krempl's  model. 

To  illustrate  this  point  an  example  problem  is  now  considered.  Sev¬ 
eral  uniaxial  bars  composed  of  mild  steel  are  subjected  to  constant  strain 
rates  at  room  temperature  with  material  properties  as  described  in  Table 
2  of  reference  [11].  Stress-strain  behavior  and  resulting  temperature 


rise  are  shown  in  Figs.  2  and  3.  These  results  were  obtained  by  integrating 
equations  (31)  and  (35)  with  a  stable  and  accurate  Euler  forward  integra¬ 
tion  scheme.  Due  to  the  rate  insensitive  nature  of  mild  steel  at  room  tem¬ 
perature,  the  predicted  results  are  identical  for  strain  rates  ranging  from 
0.001  sec  *  to  1.0  sec  *.  The  negligence  of  heat  flux  over  such  a  wide 
range  of  strain  rates  is  valid  only  under  adiabatic  conditions. 

It  is  significant  to  note  that  the  results  obtained  in  Figs.  2  and 
3  are  identical  to  those  obtained  by  Cernocky  and  Krenpl  [11].  This  is 
due  to  the  fact  that  the  assumed  internal  energy  rate  described  by  equa¬ 
tion  (14)  in  reference  [11]  can  be  obtained  in  the  uniaxial  form  by  util¬ 
izing  equations  (19),  (21)  and  (26)  in  this  paper.  Further,  energy  balance 
equation  (55)  in  reference  [11]  can  be  shown  to  be  identical  to  equation 
(26)  derived  herein  by  substituting  equation  (32)  in  this  paper.  Finally, 
it  should  be  pointed  out  that  under  non-adiabatic  conditions  neglecting 
the  heat  flux  in  the  results  obtained  herein  causes  increasing  overesti¬ 
mation  of  the  temperatures  shown  in  Fig.  3  as  the  input  strain  rate  decreases. 

Bodner  and  Partom's  model  [35,36]  may  also  be  written  in  the  uniaxial 
form  described  by  equation  (32) ,  where 


where  Dq  and  n  are  experimentally  obtained  material  constants  and 


(37) 


(38) 


where  is  an  internal  state  variable  representing  drag  stress  and  m, 

Zj,  Z^,  A,  and  r  are  experimentally  determined  material  constants.  Although 
equation  (38)  contains  stress  o ,  it  can  be  written  in  the  form  described 
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by  equation  (18)  by  direct  substitution  of  equation  (32).  Thus,  Bodner 
and  Partom’s  model  contains  two  internal  state  variables  in  the  form  de¬ 
scribed  above. 

Results  are  shown  in  Figs.  4  through  7  for  uniaxial  bars  of  IN 
100  pulled  at  various  constant  strain  rates  at  an  initial  temperature  of 
1005°K  (1350°F).  Experimental  data  were  obtained  from  Reference  [37], 
and  the  material  constants  described  above  are  tabulated  in  Reference  [38]. 
The  stress-strain  curves  shown  in  Fig.  4  are  identical  to  those  previously 
obtained  [38]. 

As  a  second  example  using  Bodner  and  Partom's  model  a  uniaxial  bar 
of  IN  100  with  material  parameters  as  described  in  Refs.  [37]  and  [38] 
is  subjected  to  the  cyclic  strain  history  shown  in  Fig.  8  and  at  initial 
temperature  1005°K  (1350°F) .  Analytic  stress-strain  behavior  is  compared 
to  experiment  in  Fig.  8  and  predicted  temperature  change  is  shown  in  Fig.  9. 

Finally,  a  uniaxial  bar  is  subjected  to  the  multicycle  test  described 
in  Fig. 10,  with  resulting  temperature  rise  shown  in  Fig.  11  It  is  observed 
that  the  model  predicts  a  mean  temperature  rise  of  approximately  3.7°K 
(6.7°F)  per  cycle.  The  linear  increase  in  mean  temperature  with  time  is 
predicted  due  to  the  cyclic  saturation  of  the  material  on  the  second  cycle, 
which  is  in  agreement  with  experimental  observations  at  elevated  temperature. 


Conclusion 

A  model  has  been  presented  herein  for  predicting  the  temperature  rise 
in  uniaxial  bars  composed  of  thermoviscoplastic  metallic  media.  The  model 
is  also  applicable  to  multiaxial  conditions,  and  this  has  been  reported  to 
some  extent  in  reference  [20],  Although  the  procedure  used  here  differs 
from  that  proposed  in  reference  [11],  it  has  been  shown  that  the  predicted 
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Fig.  5.  Predicted  Temperature  Change  for  IN100  at  1005 °K  (1350*1 
Subjected  to  Load  Histories  Shown  In  Fig.  •* 
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temperature  change  is  identical  to  results  obtained  by  Cernocky  and  Krempl 
when  their  mechanical  constitutive  equations  are  used.  However,  it  has  been 
shown  herein  that  the  introduction  of  internal  state  variables  leads  to  a 
more  general  model  which  may  be  used  with  virtually  any  thermoviscoplastic 
model  currently  used  for  metals  [24]. 

It  has  been  found  in  the  current  research  that  significant  heating  may 
occur  under  adiabatic  conditions,  expecially  during  cyclic  loading,  in  ther¬ 
moviscoplastic  metallic  media.  The  significance  of  this  heating  is  compounded 
by  the  fact  that  material  properties  often  become  extremely  sensitive  in  the 
inelastic  range  of  behavior.  This  issue  has  not  been  considered  herein,  but 
it  certainly  warrants  study  when  transient  temperature  models  become  available. 

Two  important  questions  have  not  been  answered  in  this  research:  1) 
what  effect  does  the  inclusion  of  the  heat  flux  term  have  on  the  predicted 
results;  and  2)  what,  if  anything,  does  the  present  model  have  to  do  with 
experimentally  observed  results?  The  first  question  can  only  be  addressed 
if  spacial  variation  is  admitted  in  the  field  parameters.  The  author  is  cur¬ 
rently  studying  this  question  and  hopes  to  present  results  in  a  future  paper. 
The  second  question  cannot  be  answered  at  this  time  since  it  requires  ex¬ 
tremely  sophisticated  experimentation.  Although  experimental  results  have 
been  obtained  detailing  heat  generation  in  inelastic  media,  it  is  not  pos¬ 
sible  to  compare  the  current  model  since  additional  complex  tests  must  be 
performed  in  order  to  characterize  the  thermoviscoplastic  mater?'',  parameters. 
The  author  also  hopes  to  address  this  issue  in  a  future  paper. 
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INTRODUCTION 


In  recent  years  the  analysis  of  elastic-plastic  solids  which  behave 
according  to  classical  rate-independent  incremental  plasticity  consti¬ 
tutive  models  has  become  quite  commonplace.  By  far  the  most  often  used 
method  for  nonhomogeneous  boundary  value  problem  solutions  is  the  finite 
element  method.  By  the  nature  of  the  kinematics  and  material  behavior 
this  is  a  nonlinear  field  problem,  and  a  considerable  body  of  research 
has  been  generated  dealing  with  efficient  numerical  solution  of  the 
nonlinear  boundary  value  problem.  One  efficiency  measure  adopted  solely 
to  improve  material  models  and  independent  of  large  deformation  behavior 
is  the  subincrementation  method  [1-3]. 

In  this  paper  the  subincrementation  method  will  be  briefly  reviewed 
and  an  alternative  procedure  will  be  proposed.  It  will  be  shown  that 
this  new  method  gives  both  improved  accuracy  and  efficiency  over  the  sub¬ 
incrementation  method,  and  this  contention  will  be  supported  via  several 
numerical  results. 

REVIEW  OF  THE  ELASTIC-PLASTIC  FIELD  PROBLEM 

An  elastic-plastic  medium  subjected  to  an  isothermal  loading  must 
in  general  satisfy  the  following  conditions  at  all  points  in  its  interior 
V  and  on  its  surface  B:  (1)  conservation  of  linear  and  angular  momentum, 
(2)  conservation  of  mass;  (3)  strain-displacement  relations;  (4)  the 
first  and  second  laws  of  thermodynamics;  and  (5)  stress-strain  relations 
(constitution).  Of  course,  there  are  additional  constitutive  relations 
[4],  but  these  need  not  be  considered  in  order  to  characterize  the  mech¬ 
anical  response  when  internal  heat  generation  is  negligible.  It  will 


be  assumed  that  the  above  restriction  holds  for  the  media  considered 
herein,  although  the  method  to  be  considered  here  applies  equally  to  tran 
sient  temperature  phenomena. 

Condition  (5)  is  the  main  topic  of  the  discussion  herein.  In  this 
paper  we  consider  only  elastic -plastic  media:  that  class  of  materials 
for  which  the  stress  (or  strain)  tensor  can  be  considered  to  be  a  time 
independent  functional  of  the  strain  (or  stress)  tensor,  that  is,  stress 
is  dependent  on  the  entire  history  of  strain  but  independent  of  the  time 
scale.  It  has  been  shown  that  this  functional  form  can  often  be  written 
in  an  equivalent  differential  equation  form  [5].  One  common  strain  for¬ 
mulation  is  of  the  general  type  [6] 


da.  .  -  C.  .de. 
ij  ljkl  kl 


where  for  small  deformation  o^_.  and  represent  infinitesimal  stress 
and  strain  tensors,  respectively,  and  called  the  effective  mod¬ 

ulus  tensor,  is  the  repository  for  history  dependence  via  its  dependence 

_P 

on  the  equivalent  uniaxial  plastic  strain,  e  ,  which  is  a  metric  in  the 
space  of  plastic  strain  defined  by 


eP(t) 


eP(t)  = 


eP  (t) 

n 


2  ,  p  ,  P 
r  de. .de.  . 

3  ij  ij 


where  t  represents  time  and 


dCij  =  dCij  '  Dijkldakl 


and  D  *s  the  elastic  modulus  tensor,  assumed  to  be  independent  of 
deformation.  In  a  uniaxial  sense  the  effective  modulus  tensor  C,.,  , 


may  be  thought  of  as  a  secant  modulus. 


It  should  be  noted  that  for  finite  deformation  equations  (1)  may 
still  be  applicable  if  and  are  replaced  by  frame  indifferent 
quantities  consistent  with  hypoelasticity  [7].  Note  also  that  under 
certain  simplifying  assumptions  equations  (1)  reduce  to  the  classical 
Prandtl- Reuss  equations  [8,9]. 

It  has  been  shown  that  equations  (1)  are  consistent  with  the  first 
and  second  laws  of  thermodynamics  {10-12].  Therefore,  since  the  energy 
balance  is  trivially  satisfied  in  an  isothermal  domain  with  negligible 
internal  heat  generation  and  conservation  of  mass  is  satisfied  if  the 
density  is  timewise  constant  during  infinitesimal  deformation,  we  need 
consider  only  conditions  (1),  (3)  and  (5)  here. 

In  the  finite  element  method  conditions  (1)  and  (3)  are  usually 
satisfied  via  an  incremental  variational  principle  integrated  over  the 
domain  of  interest.  The  discretization  process  then  entails  reducing 
the  volume  integral  for  the  variational  principle  to  some  sub-domain 
aptly  called  a  finite  element.  The  integration  over  the  volume  of  each 
element  is  usually  sufficiently  difficult  to  require  numerical  integra¬ 
tion,  and  for  this  purpose  a  quadrature  procedure  is  generally  employed. 
Thus,  the  integration  is  reduced  to  evaluation  of  the  integrand  at  a 
finite  number  of  integration  points  within  an  element. 

The  volume  integration  of  each  element  yields  a  set  of  matrix  equa¬ 
tions  which  are  assembled  into  a  global  set  of  matrix  equations.  These 
equations  are  nonlinear  even  during  infinitesimal  deformations  due  to 
the  nonlinearity  of  equations  (1).  Therefore,  an  iterative  technique 
is  generally  used  to  solve  the  global  equations  on  each  load  step,  and 
the  constitutive  equations  (1)  must  be  solved  several  times  for  all 


integration  points  on  each  load  step  until  convergence  occurs.* 


SOLUTION  OF  THE  CONSTITUTIVE  EQUATIONS 


It  has  become  common  practice  in  the  literature  to  incrementalize 
equations  (1)  by  simply  replacing  differentials  da^ ^  and  d£„  with  in¬ 
crements  Aa_  and  ,  respectively.  However,  since  the  integrand 
depends  on  the  strain  tensor  during  the  load  increment,  an  important 
task  becomes  the  integration  of  equations  (1)  over  some  input  increment 
in  the  strain  tensor,  viz.: 


o. . (e  . (t+At))  e. , Ct+At) 


A  a.  .  ■ 


J  d°ij  ~  J  Cijkl  (£  )dCkl 

°u(eij(t)) 


where  t  represents  the  time  at  the  start  of  a  load  increment,  and  t+At 
is  the  time  at  the  end  of  a  load  increment. 

Equations  (4)  definitely  present  a  uniqueness  problem  since  the 
strain  tensor  may  be  cycled  during  the  time  increment  At.  In  order  to 
avoid  this  difficulty  a  sufficient  condition  may  be  adopted  which  is 
not  unlike  the  condition  required  to  obtain  the  Mises-Hencky  deforma¬ 
tion  theory  [15]  from  the  Prandtl-Reuss  equations.  It  is  assumed  that 
during  the  time  increment  At  all  components  of  the  strain  tensor  increase 
monotonically  via  the  relation 

— p 

dekl  «=  K^de  ,  =  constants  ,  (5) 

— P 

where  de  must  be  a  monotonically  increasing  function  of  strain  during 
plastic  loading  over  the  time  increment  At.  Substitution  of  equations 
(5)  into  (4)  yields 


*For  a  more  complete  discussion  of  the  finite  element  method  applied  to 
elastic-plastic  media,  see  references  [13]  and  [14]. 
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which  is  obviously  a  unique  relation  due  to  the  monotinicity  of  the 
plastic  strain  increment  during  the  time  increment  At,  Now  define  the 
following  fourth  order  tensor  which  is  constant  over  the  time  step  At: 


C' 

ijkl 
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C^(t+At) 
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Then  substitution  of  definitions  (7)  into  equation  (6)  gives 


4oii  “ 
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(8) 


which  is  the  exact  incremental  relation  which  should  be  used  with  the 
incremental  variational  principle.  It  should  be  noted  that  equations 


(1)  and  (8)  are  by  no  means  equivalent  since  can  seen  from 


definitions  (7)  to  represent  an  average  effective  modulus  tensor  during 


the  time  increment  At.  Unfortunately,  equations  (7)  cannot  be  integrated 

HP, 


precisely  because  the  upper  limit  of  integration  c  (t+At)  cannot  be  de¬ 
termined  until  equations  (8)  have  been  evaluated. 

It  will  be  recalled  that  the  equivalent  uniaxial  plastic  strain 
can  be  shown  to  be  equal  to  the  axial  plastic  strain  when  a  bar  is  pulled 
uniaxially  [16].  Now  define  the  equivalent  uniaxial  stress 


0  = 


*■[? 
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which  can  also  be  seen  to  be  equivalent  to  the  axial  stress  when  a  bar 
is  pulled  uniaxially  [16].  Thus,  the  information  required  to  charac¬ 
terize  equations  (6)  is  obtainable  from  a  single  raonotonically  increasing 


-  6  - 


equivalent  uniaxial  plastic  strain  diagram  as  shown  in  Fig.  1.  Note 
that  the  curve  shown  represents  nothing  more  than  a  uniaxial  stress- 
strain  diagram  with  abscissa  transformed  via  definition  (3).  In  addi¬ 
tion,  for  combined  isotropic-kinematic  hardening  the  ordinate  in  Fig. 

1  should  be  transformed  as  well  [17]. 

It  is  apparent  from  Fig.  1  that  for  a  continuously  work  hardening 

material  the  relation  between  the  equivalent  uniaxial  plastic  strain 
— P 

e  and  the  slope  of  the  uniaxial  stress-plastic  diagram  is  unique,  that 


■'(f) 


where  F  is  a  bijective  mapping.  Therefore,  the  effective  modulus  tensor 


may  be  written  alternatively  as 


(§)  ■ 


■  Vi<eP) 


Thus,  because  the  effective  modulus  tensor  C...,  is  a  nonlinear  function 

ljkl 

_  _ p 

of  the  plastic  secant  modulus  dG/dE  ,  integration  of  equations  (7)  is 
not  a  trivial  task.  In  order  to  avoid  this  numerical  complexity  it  is 
not  uncommon  to  simply  approximate  the  effective  modulus  tensor  by 

cijki^P)  !  cijki(7P(t))  •  (12> 

thus  reducing  integration  of  equations  (6)  to 


ioij  *  cijki<cP(t»\i4cP  ■  cijki^P<t»Acki  •  <13: 

which  is  Euler's  method  of  forward  integration.  Obviously,  since  this 
is  nothing  more  than  a  simple  first  order  Taylor  series  expansion  its 
accuracy  will  depend  on  the  relative  nonlinearity  or  curvature  of  the 
uniaxial  stress-plastic  strain  curve  during  a  given  load  increment. 
This  condition  is  illustrated  for  the  uniaxial  case  in  Fig.  2. 


—APPROXIMATE 


Fig.  2. 
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In  many  computer  codes  the  uniaxial  stress-strain  diagram  is  input 
in  piecewise  linear  fashion  as  shown  in  Fig.  3.  While  it  is  not  clear 
that  this  piecewise  linearization  is  motivated  by  anything  beyond  sim¬ 
plification  of  input  data,  it  has  the  added  benefit  that  it  helps  alle¬ 
viate  the  numerical  integration  problem  noted  above  and  described  in 

— p 

Fig.  4.  In  fact,  so  long  as  the  plastic  strain  increment  A £  does  not 
subtend  a  slope  discontinuity  during  a  load  step  equations  (11)  indi¬ 
cate  that  approximation  (12)  will  reproduce  the  piecewise  linear  curve 
exactly.  Therefore,  the  accuracy  of  equations  (13)  is  limited  only  by 
the  accuracy  with  which  one  can  reproduce  the  exact  curve  of  stress  versus 
plastic  strain  with  a  piecewise  linear  curve.  Mathematically,  the  slope 

continuity  condition  is  satisfied  if  one  can  find  values  of  the  equiva- 

— P 

lent  plastic  strain  at  slope  discontinuities  E^,  as  shown  in  Fig.  3, 
such  that  for  the  current  load  increment 


(14) 


for  all  equivalent  plastic  strains  in  the  range 

7p(t)  _<  rp  _<  rp(t+At)  .  (15) 

However,  condition  (14)  cannot  be  a  priori  guaranteed  in  practice  be¬ 
cause  in  a  non-homo geneous  boundary  value  problem  the  equivalent  uniaxial 
plastic  strain  varies  spatially.  Whereas  one  integration  point  may 
undergo  a  very  small  or  even  zero  (elastic)  plastic  strain  increment 
during  a  specified  increment  in  surface  tractions,  another  point  under 
high  stress  concentration  may  undergo  a  plastic  strain  increment  which 
subtends  one  or  more  discontinuities  in  the  piecewise  linear  equivalent 
uniaxial  stress  versus  equivalent  uniaxial  plastic  strain  diagram.  Thus, 
as  can  be  seen  from  definitions  (2)  and  (3)  the  equivalent  uniaxial 


Fig.  3.  Piecewise  Linearization  of  the  Equivalent 
Uniaxial  Stress  Versus  Equivalent  Uniaxial  Plastic 
Strain  Diagram 
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Fig.  4.  Estimated  Equivalent  Uniaxial  Stress 
Increment  Using  Subincrementation 


plastic  strain  increment  can  be  determined  only  after  the  stress  incre¬ 
ment  tensor  has  been  found. 


REVIEW  OF  THE  SUBINCREMENTATION  METHOD 


In  order  to  improve  the  accuracy  of  approximation  (13)  subincremen¬ 
tation  has  been  proposed  [1,2].  In  this  method  Stricklin,  et  al.  define 
the  equivalent  uniaxial  strain  increment 


d£  -  V  3  dEijdEij  *  (16) 

This  quantity  is  evaluated  over  a  specified  load  step  and  is  then 
compared  to  an  input  parameter  called  the  allowable  total  strain  incre¬ 
ment  (de._)  as  follows 
AL 


where  M  is  rounded  off  to  the  nearest  integer  greater  than  zero.  Equa¬ 
tions  (6)  are  then  evaluated  M  times  for  the  strain  subincrement 

S  ^Cii 

5  -r1  •  <18> 

and  on  each  subincrement  the  effective  modulus  tensor  is  updated  to 
reflect  the  current  equivalent  uniaxial  plastic  strain.  Based  on  numer¬ 
ical  evidence,  Stricklin  suggests  that  de^  should  be  no  greater  than 
0.0005  in. /in. ,  although  our  experience  indicates  that  values  as  small 
as  .00005  in. /in.  are  sometimes  required  to  maintain  accuracy  of  solu¬ 
tion. 


In  order  to  illustrate  the  effect  of  subincrementation  let  us  ex¬ 


amine  a  single  example.  Suppose  we  consider  a  bar  subjected  to  a  grad¬ 
ually  increasing  homogeneous  uniaxial  stress  state.  Because  conditions 


(1)  and  (3)  are  satisfied  trivially  we  need  only  consider  approximate 
constitutive  equations  (13).  Since  the  input  material  properties  will 
be  described  via  a  piecewise  linear  equivalent  uniaxial  stress  versus 
equivalent  uniaxial  plastic  strain  diagram  as  shown  in  Fig.  3,  and  be¬ 
cause  this  boundary  value  problem  is  equivalent  to  the  experiment  which 
produced  the  material  input  data,  an  exact  analysis  using  equations 
(13)  should  reproduce  Fig.  3  precisely.  In  fact,  using  subincrementa¬ 
tion  will  yield  the  results  shown  in  Fig.  4  when  a  single  slope  discon¬ 
tinuity  is  encountered  in  a  given  load  step.  It  can  be  seen  from  the 
figure  that  the  total  error  is  incurred  during  the  plastic  strain  sub¬ 
increment  subtending  the  slope  discontinuity.  The  effect  then  of  sub¬ 
incrementation  is  simply  to  improve  the  approximate  integration  of  equa¬ 
tions  (6). 

A  PROPOSED  MODIFICATION 

It  is  apparent  from  the  above  discussion  that  subincrementation 
will  often  require  multiple  evaluations  of  equations  (13)  for  each  inte¬ 
gration  point.  Since  these  equations  must  be  evaluated  at  each  inte¬ 
gration  point  in  the  body  and  often  several  times  for  each  load  step 
in  order  to  obtain  equilibrium  convergence,  considerable  computational 
time  can  be  spent  in  this  process.  Detailed  herein  is  a  numerical  pro¬ 
cedure  for  integrating  equations  (7)  which  is  both  more  accurate  and 
more  computationally  efficient  than  subincrementation.  We  call  this 
method  the  zeta  method. 

The  method  proposed  here  is  a  simple  extension  of  a  procedure  uti¬ 
lized  by  Krieg  and  Duffey  (18]  for  the  transition  step  from  elastic  to 
elastic-plastic  material  behavior.  The  primary  extension  is  that  each 


subsequent  slope  discontinuity  in  the  piecewise  linear  equivalent  uni¬ 


stress  versus  equivalent  uniaxial  plastic  strain  diagram  is  treated  ex¬ 


actly  like  a  subsequent  yield  surface.  Since  equations  (13)  are  exact 


under  conditions  (1A)  and  (15),  no  subincrementation  will  be  required 


to  obtain  precise  results  between  slope  discontinuities. 


In  order  to  see  how  the  zeta  method  works,  consider  a  material  point 


which  is  in  a  post-yielded  state  at  time  t  and  with  equivalent  uniaxial 


— p 

plast  strain  £  (t) ,  as  shown  in  Fig.  5. 


According  to  Krieg  and  Duffey  [18],  the  value  of  the  stress  tensor 


at  the  material  point  necessary  to  bring  the  equivalent  uniaxial  stress 


state  to  the  i+lth  slope  discontinuity  is  defined  by 


♦  C  AS 


where  Ao^  is  the  increment  predicted  by  equations  (13)  using  the  input 


total  strain  increment  Ae^  and  £  is  a  scalar  factor  to  be  determined. 


In  order  to  determine  zeta  definitions  (19)  are  substituted  into 


the  yield  criterion  used  in  the  model.  For  example,  if  von  Mises' 


yield  criterion  is  utilized,  equations  (19)  will  result  in 


(owHtAow) 


(0t„,+LAo„T,) 


Solving  the  above  equation  for  zeta  will  result  in 


-B  4-  y  (B  -  AAC) 


where 


A  -  Ao^Ao^  » 


B  -  20, j  iOy  - 


-  15  - 


v' 


WWW 
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and 


c  *  °yaij  -  (0i+i>2  •  (24> 

Utilizing  equations  (21)  through  (24)  the  value  of  zeta  may  be 
obtained.  If  zeta  is  greater  than  or  equal  to  unity,  the  input  total 
strain  increment  will  not  take  the  equivalent  uniaxial  stress  beyond 
the  next  slope  discontinuity  and  the  results  of  equations  (13)  may 
be  considered  correct.  In  other  words 

£  >  1  =>  A o. .  »  Ao. .  .  (23) 

-  iJ 

If,  on  the  other  hand,  zeta  is  less  than  unity,  the  predicted  stress 
increment  tensor  is  incorrect  and  equations  (13)  must  be  modified. 

This  is  accomplished  by  first  constructing  the  input  strain  increment 
necessary  to  bring  the  equivalent  uniaxial  stress  to  the  slope  discon¬ 


tinuity: 


.  i+1 


where  Ae  is  the  input  strain  increment.  The  values  of  Ae . .  are 
ij  ij 

then  substituted  into  equations  (13)  to  produce 

Ao1+1  •=  C.  ..  .  (IP(t))Ae*+1  ,  (27) 

1J  ljkl  '  ij 

the  remaining  portion  of  the  stress  increment  tensor  is  calculated 
by  first  determining  the  remainder  of  the  input  strain  increment 

AeJ.  =  (1  -  OAe.j  ,  (28) 

noting  from  definitions  (26)  and  (28)  that 

Si  ■  <+1  +  ^  •  <29) 

R 

The  remainder  of  the  strain  increment  tensor  Ae^  then  substituted 
into  equations  (13)  to  give  the  remainder  of  the  stress  increment  tensor 


17  - 


I'.TT.TV’M  VV  1%  W 


.  R  _  P  ..  R 

A°ij  *  CijklCci+l)A  ij 


Thus,  the  total  stress  increment  tensor  is  given  by 


Ao  =  Aoi+1  +  Ao?. 
ij  ij  ij 


It  is  easily  verified  that  the  above  procedure  will  result  in  an  equiv¬ 


alent  uniaxial  stress  and  equivalent  uniaxial  plastic  strain  [utilizing 


definitions  (2)  and  (3)]  which  lie  on  the  equivalent  uniaxial  stress 


versus  equivalent  uniaxial  plastic  strain  diagram.  It  should  also  be 


pointed  out  that  although  the  actual  yield  surface  is  updated  through¬ 


out  plastic  loading,  the  equivalent  uniaxial  stresses  corresponding  to 


slope  discontinuities  should  at  no  time  be  altered. 


Although  the  above  procedure  has  been  discussed  here  only  in  the 


context  of  isotropic  hardening,  it  is  also  applicable  to  more  complex 


yield  criteria  and  work  hardening  rules  [19,20], 


COMPUTER  CODE  FLOWCHART 


The  following  chart  outlines  in  abbreviated  form  the  application 


of  the  <5  method  for  a  given  increment  in  the  total  strain  tensor  Ac,. 

ij 


and  equivalent  uniaxial  plastic  strain  at  the  start  of  the  step  er(t). 


a)  Set  Ao . .  *  0 , 
ij 


b)  Evaluate  C^ ^  (c  (t)). 


c)  Obtain  Ao^  using  equations  (13). 


d)  Determine  from  Fig.  5. 


e)  Calculate  A,  B,  and  C  using  equations  (22)  through  (24). 


f)  Determine  zeta  using  equation  (21). 


g)  If  (  >  1  go  to  step  q). 


-v-  (v \v..' 
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h)  If  £  <  1  evaluate  Ae^  using  equations  (26), 


i)  Determine  Aa^^  using  equations  (27). 

j)  Set  A a. 


a  .  a  i+l 
'a  M  4oij  +  ioij  • 


k)  Calculate  Ae„  using  equations  (28), 

l)  Set  Ae . .  *  Aef.. 

iJ  iJ 

P 

m)  Calculate  Ae_  using  equations  (3). 

— p 

n)  Determine  Ae  using  equation  (2). 

o)  Set  eF(t)  *  eP(t)  +  AeP. 

p)  Go  to  step  b). 

1)  Set  Aa.  .  «  Ao.  .  +  A oJ  .. 

ij  iJ  ij 


DISCUSSION  OF  RESULTS 


In  this  section  the  results  generated  using  the  zeta  method  as 
well  as  the  subincrementation  method  for  solving  the  constitutive  equa¬ 
tions  of  classical  plasticity  are  presented.  Both  solution  techniques 
have  been  incorporated  into  a  finite  element  computer  program  which 
uses  constant  strain  triangular  elements.  The  formulations  have  been 
cast  into  a  2-dimensional  plane  stress  format. 

In  order  to  compare  the  efficiency  of  the  two  different  methods 
their  respective  solution  times  will  be  compared.  This  was  accomp¬ 
lished  by  using  the  built  in  timer  (clock)  subroutine  used  in  the 
Fortran-H  Extended  language  available  on  the  AMDAHL  470  V6  located 
on  the  Texas  A&M  University  campus.  Two  times  will  be  given  in  the 
analysis:  1)  time  spent  using  the  constitutive  package,  and  2)  total 

time  spent  in  solving  the  specified  boundary  value  problem. 

Two  boundary  value  problems  have  been  selected  for  comparing  the 


two  constitutive  packages:  a  highly  yielded  uniaxial  bar  subjected 


to  uniaxial  tension  only  and  a  thick  walled  pressure  vessel  subjected 
to  internal  pressure  sufficient  to  yield  a  broad  band  of  the  vessel. 
Both  specimens  are  assuned  to  be  made  of  5086-H34  aluninun  with  piece- 
wise  linearized  room  temperature  properties  shown  in  Figure  6.  The 
finite  element  mesh  used  for  the  uniaxial  bar  as  well  as  the  load 
input  diagram  are  shown  in  Figure  7.  Only  two  elements  are  necessary 
to  represent  the  bar  because  the  boundary  value  problem  is  homogeneous. 
However,  if  the  problem  was  inhomogeneous  then  mesh  refinement  would 
be  necessary  to  increase  the  accuracy  of  the  solution.  It  should 
nevertheless  be  pointed  out  here  that  the  number  of  elements  used  in 
the  mesh  is  directly  proportional  to  the  computational  time  required 
in  the  constitutive  package.  Another  factor  influencing  the  computa¬ 
tion  time  is  the  non-linearity  of  the  given  stress  strain  curve.  In¬ 
creasing  the  nimber  of  piecewise  linearities  in  the  idealized  uniaxial 
stress  strain  curve  will  increase  the  computational  time  required  by 
the  zeta  method.  Although  this  increase  in  piecewise  linearizations 
will  not  greatly  affect  the  computational  time  required  by  the  subin¬ 
crementation,  it  will  have  an  adverse  affect  on  the  accuracy  of  this 
method. 

The  results  of  the  uniaxial  bar  test  are  shown  in  Fig.  8  and  the 
comparative  solution  times  are  given  in  Tables  1  and  2.  Fig.  8  shows 
the  output  axial  displacement  versus  time  for  the  zeta  method  as  well 
as  the  subincrementation  with  various  allowable  errors  in  equivalent 
uniaxial  strain  shewn  in  parentheses.  Table  1  shows  a  comparison  of 
solution  times  for  each  load  step,  while  Table  2  gives  a  more  detailed 
comparison  of  solution  times  for  the  final  time  step.  Several  dif¬ 
ferent  cases  were  run  using  the  subincrementation  code  in  order  to 


UNIAXIAL  STRESS  STRAIN  CURVE 


UNIAXIAL  BAR  TENSION  TEST 

5D86-H34  ALUMINUM 
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TABLE  2: 


COMPARISON  OF  SOLUTION 

TIMES  FOR 

THE 

LAST  LOAD  STEP 

ON  THE 

UNIAXIAL  BAR 

SOLUTION  TIMES 

(sec) 

TIME 

ZETA 

MET  HDD 

SUBINC 

(0.005) 

SUBINC 

(0.0005) 

SUBINC 

(0.0001) 

SUBINC 

(0.00005) 

CONSTITUTIVE 

0.063492 

0.0027 

0.0555 

0. 18764 

0.37684 

SUBINC/.  063492 

.043 

.716 

2.955 

5.935 

TOTAL  TIME  (SEC) 

0.0968 

0.00798 

0.07761 

0.22315 

0.40994 

SUBINC/.  0968 

.0824 

.802 

2.305 

4.  235 

illustrate  the  difference  in  solutions  and  comparative  times  for  dif¬ 
ferent  values  of  the  allowable  uniaxial  stress  increment  de.T  used 

AL 

in  determining  the  number  of  subincrements.  In  general,  increasing 
the  required  accuracy  of  the  subincrementation  method  also  increases 
the  computational  time. 

The  results  of  this  test  show  that  in  order  to  obtain  solution 

accuracy  by  subincrementation  which  is  comparable  to  the  zeta  method, 

the  allowable  strain  increment  must  be  of  the  order  of  0.00005  in/in. 

In  fact,  the  larger  the  allowable  strain  increment  the  less  solution 

time  required,  and  in  fact  for  dt.,  *  0.005  and  0.0005  the  solution 

AL 

time  was  less  for  the  subincrementation  method.  However,  the  resulting 
accuracy  was  very  poor.  Table  1  indicates  that  subincrementation 
(e^  *  0.00005  in. /in. )  requires  3.171  times  as  much  computer  time  as 
the  zeta  method  for  the  uniaxial  bar  problem.  Although  this  is  a 
rather  large  difference  in  relative  times,  since  only  a  two  element 
problem  has  been  run,  the  difference  in  actual  cost  is  small.  However, 
on  an  extremely  large  scale  problem  obvious  savings  would  result. 

The  finite  element  mesh  used  for  the  thickwall  pressure  vessel 
is  shown  in  Fig.  9  and  the  load  input  diagram  is  shown  in  Fig.  10. 

The  results  of  this  test  are  shown  in  Fig.  1 1  as  well  as  Tables  3  and 
A.  These  results  are  for  the  final  pressure  of  P  *  30,000  PSI.  It 
should  be  noted  that  if  one  applies  increasing  pressure  to  the  speci¬ 
men  then  more  of  the  elements  will  yield  in  the  outer  regions  of  the 
thickwalled  pressure  vessel,  resulting  in  a  higher  solution  time  be¬ 
cause  more  time  is  spent  in  the  constitutive  package  and  more  time 
in  iterating  on  the  correct  nonlinear  solution. 


PRESSURE  (xlOOO  lbs) 


THICK  WALLED  PRESSURE  VESSEL 


RADIAL  LOCATION  IN 


TABLE  3:  RADIAL  DISPLACEMENTS  AND  SOLUTION 

TIMES  FOR  THICK  WALLED  PRESSURE  VESSEL 


RADIAL  RADIAL  DISPLACEMENT  (IN) 

LOCATION 


ZETA 

SUBINC 

SUBINC 

SUBINC 

SUBINC 

METHOD 

(0.001) 

(0.0005) 

(0.0001) 

(0.00005) 

2.0 

0.01324344 

0.01285087 

0.01289475 

0.01292329 

0.0129211 

2.2222 

0.01167517 

0.01133769 

0.01137959 

0.01141027 

0.01140665 

2.4444 

0.0105329 

0.01026293 

0.01029393 

0.01032073 

0.01031638 

2. 6666 

0.00971329 

0.00947785 

0.00960516 

0.00953320 

0.00952741 

3.0 

0.00883162 

0.00862648 

0.00866162 

0.00868949 

0.00868414 

3.3333 

0.00821108 

0.00802864 

0.00806262 

0.00809172 

0.00808688 

4.0 

0.00728626 

0.00714737 

0.00717933 

0.00720508 

0.00720022 

CONSTITUTIVE 

PACKAGE  TIME  (SEC) 

0.691703 

0.6132077 

0.6854118 

1.056744 

1.503111 

SUBINC/. 691703 

.886 

.991 

1.527 

2.173 

TOTAL  TIME  (SEC) 

3.1546827 

3.0397627 

3.1661748 

3.61964 

4.0580 

SUBINC/3. 1546827 

.  964 

1.003 

1.147 

1.286 
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TABLE  4: 


COMPARISON  OF  SOLUTION  TIMES  FOR 
THE  LAST  LOAD  STEP  ON  THE 
THICK  WALLED  PRESSURE  VESSEL 


SOLUTION  TIMES  (sec) 


TIMES 

ZETA 

MET  TOD 

SUBINC 

(0.001) 

SUBINC 

(0.0005) 

SUBINC 

(0.0001) 

SUBINC 

(0.00005) 

CONSTITUTIVE 

PACKAGE  TIME  (SEC) 

0.476943 

0.398788 

0.468156 

0.838656 

1.278211 

SUBINC/.  476943 

.830 

.  982 

1.758 

2.680 

TOTAL  TIME  (SEC) 

1.662752 

1.565564 

1.665976 

2.143309 

2.579693 

SUBINC/1. 662752 

.941 

1.002 

1.289 

1.551 

The  results  of  the  test  are  basically  the  same  as  those  for  the 
uniaxial  bar.  From  Table  3  it  can  be  seen  that  for  comparable  solu¬ 
tion  accuracy  the  subincrementation  method  takes  1.527  to  2.173  times 

— p 

more  constitutive  time  (depending  on  c  )  than  the  zeta  method  and 

AL 

1.147  to  1.286  times  greater  total  time.  Closer  examination  shows 

— p 

that  even  for  the  smallest  eAT  =  0.00005  the  solutions  still  differ 

AL 

from  the  zeta  method  and  that  there  is  no  noticeable  difference  be- 

— p  — p 

tween  solutions  for  de.,  «  0.0001  and  de4T  =  0.00005.  In  fact,  the 

AL  AL 

— P 

results  tend  to  be  less  accurate  for  de.,  *  0.00005.  This  can  be  at- 

AL 

tributed  to  numerical  roundoff  error  because  the  increments  in  the 
strain  are  so  small  that  further  improved  accuracy  is  not  possible. 

By  constrast,  there  is  no  nunerical  roundoff  error  apparent  in  the 
zeta  method. 

CONCLUSION 

The  objective  of  this  resear c  has  been  to  produce  an  alterna¬ 
tive  to  subincrementation  which  results  in  a  substantial  improvement 
in  computational  efficiency.  This  new  method  has  been  shown  by  ex¬ 
ample  to  give  not  only  improved  efficiency,  but  also  slightly  greater 
accuracy  of  results.  The  following  general  conclusions  can  be  made: 

1)  in  order  to  produce  results  by  the  subincrementation  method 
which  are  comparable  in  accuracy  to  the  zeta  method,  signi¬ 
ficantly  greater  computation  time  is  required  by  the  former 
method; 

2)  increasing  required  accuracy  in  the  allowable  equivalent  uni¬ 
axial  strain  increment  (Ac.,)  can  lead  to  roundoff  error  when 

AL 

subincrementation  is  utilized; 


V-l-V'v 


V" V*  v.*  o  «_•  * 


3)  subincrementation  necessarily  produces  errors  in  predicted 
stresses  whenever  a  slope  discontinuity  is  subtended  in  the 
uniaxial  stress-strain  curve; 

4)  the  zeta  method  follows  the  uniaxial  stress-strain  curve  ex¬ 
actly; 

5)  both  subincrementation  and  the  zeta  method  approximate  the 
load  path  to  be  radial  during  each  load  increment, 

6)  piecewise  linearization,  although  merely  a  numerical  conven¬ 
ience  in  subincrementation,  is  necessary  in  order  to  utilize 
the  zeta  method; 

7)  the  zeta  method  can  be  used  with  cyclic  hardening  models  of 
plasticity;  and 

8)  the  zeta  method  may  not  be  appropriate  for  use  in  rate  depen¬ 
dent  viscoplastic  media. 
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ABSTRACT 


The  thermomechanical  response  of  a  uniaxial  bar  with  thermoviscoplastic 
constitution  is  predicted  herein  using  the  finite  element  method.  After  a 
brief  review  of  the  governing  field  equations,  variational  principles  are  con¬ 
structed  for  the  one  dimensional  conservation  of  momentum  and  energy  equations. 
These  equations  are  coupled  in  that  the  temperature  field  affects  the  displace¬ 
ments  and  vice  versa. 

Due  to  the  differing  physical  nature  of  the  temperature  and  displacements, 
first  order  and  second  order  elements  are  utilized  for  these  variables,  respec¬ 
tively.  The  resulting  semi-discretized  equations  are  then  discretized  in  time 
using  finite  differencing.  This  is  accomplished  by  Euler's  method,  which  is 
utilized  due  to  the  stiff  nature  of  the  constitutive  equations. 

The  model  is  utilized  in  conjunction  with  stress-strain  relations  devel¬ 
oped  by  Bodner  and  Partom  to  predict  the  axial  temperature  field  in  a  bar  sub¬ 
jected  to  cyclic  mechanical  end  displacements  and  temperature  boundary  condi¬ 
tions.  It  is  found  that  spacial  and  time  variation  of  the  temperature  field 
is  significantly  affected  by  the  boundary  conditions. 

TABLE  OF  SYMBOLS 

t  -  time 

P  -  axial  internal  resultant  force 

-  axial  externally  applied  force  per  unit  length 
x  -  axial  coordinate  dimension 
O  -  axial  stress  component 
A  -  cross-sectional  area 

Tx  -  end  traction  in  units  of  force  per  unit  area 
s  -  surface  area 
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Sc  -  area  of  the  longitudinal  surface  of  the  bar 
£  -  axial  strain  component 

u  -  axial  displacement  component 

-  internal  state  variable  representing  axial  inelastic  strain 
E  -  Young's  modulus  in  the  axial  coordinate  direction 

a  -  coefficient  of  thermal  expansion  in  the  axial  coordinate  direction 
T  -  temperature 

T  -  reference  temperature  at  which  no  deformation  is  observed  at  zero  load 
R 

-  internal  state  variable  representing  drag  stress 
q  -  axial  component  of  heat  flux 

k  -  coefficient  of  axial  thermal  conductivity 
Cv  -  specific  heat  at  constant  elastic  strain 
p  -  mass  density 

r  -  internal  heat  source  per  unit  mass 
L  -  length  of  the  bar 

INTRODUCTION 

It  is  well  known  that  mechanical  and  thermodynamic  coupling  are  signif¬ 
icant  in  metallic  solids  [1-11].  The  author  has  recently  developed  a  model 
capable  of  predicting  this  coupling  effect  in  thermoviscoplastic  metals  [12]. 
In  the  previous  paper  a  cyclic  strain  control  loading  on  a  sample  of  INI 00  at 
1005°K  (1350°F)  was  used  to  predict  a  temperature  rise  of  approximately  3.7°K 
per  cycle  when  the  strain  amplitude  was  2 %  and  the  specimen  was  adiabatically 
insulated. 

The  focus  of  the  current  research  is  to  consider  the  effect  of  thermal 
boundary  conditions  on  this  same  process.  The  introduction  of  these 


conditions  causes  the  strain  and  temperature  fields  to  be  inhomogeneous  even 
though  the  stress  field  is  homogeneous  if  the  bar  is  prismatic.  This  spacial 
variation  in  the  field  variables  causes  the  process  to  be  difficult  to  model 
because  the  thermomechanical  constitutive  equations  are  highly  nonlinear  stiff 
differential  equations.  In  this  paper  the  finite  element  method  is  utilized 
to  spatially  discretize  the  dependent  variables  displacement  and  temperature, 
and  the  finite  difference  method  is  employed  for  timewise  discretization. 

This  process  results  in  a  set  of  highly  nonlinear  algebraic  equations. 

Since  the  thrust  of  this  research  is  to  obtain  accurate  results  without 
regard  to  numerical  efficiency,  the  results  are  obtained  via  the  relatively 
inefficient  but  accurate  method  of  simply  utilizing  successively  smaller  time 
steps  along  with  refined  spatial  mesh  to  obtain  a  convergent  and  therefore 
accurate  solution  for  the  temperature  and  displacement  fields  both  spatially 
and  as  a  function  of  time  for  a  cyclically  imposed  end  displacement. 

The  physical  interest  in  the  problem  is  to  determine  the  effect  of 
temperature  boundary  conditions  on  the  predicted  temperature  rise  in  a  bar  sub¬ 
jected  to  cyclic  mechanical  loading.  It  is  found  from  the  analysis  that  the 
introduction  of  these  nonadiabatic  boundary  conditions  causes  significant  axial 
temperature  gradients.  Since  nonadiabatic  conditions  cannct  be  avoided  in 
experimental  research,  it  is  concluded  that  experimental  tests  of  this  type 
should  be  viewed  with  caution  when  their  purpose  is  to  construct  constitutive 
relations. 

PROBLEM  SOLUTION 


Field  Problem  Description 


The  following  field  equations  are  given: 


a)  equilibrium  [13], 


m 


where  the  axial  resultant  P  is  defined  by 
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b)  strain-displacement  relation 


c)  thermomechanical  constitution. 
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where  z  is  the  total  number  of  internal  state  variables;  and 
d)  conservation  of  energy 


and  (6) 
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The  conservation  of  mass  is  satisfied  trivially  and  the  second  law  of  thermo¬ 
dynamics  has  been  previously  shown  to  be  satisfied  by  the  above  equations  [14- 
16].  It  should  be  noted  that  equilibrium  equation  (1)  satisfies  equilibrium 
in  the  axial  coordinate  direction  only  in  an  average  sense  over  the  cross-section. 


I 


4 


The  above  6+Z  equations  (excluding  definition  (3))  define  a  nonlinear 
initial-boundary  value  problem  (together  with  appropriate  thermal  and  mech¬ 


anical  initial  and  boundary  conditions)  in  which  the  following  dependent  vari¬ 
ables  are  sought  as  functions  of  x  and  t:  0,  e,  u,  q,  T,  P,  and  ot^. 

For  convenience  the  domain  is  defined  to  be  of  length  L,  so  that  boundary 
and  initial  conditions  are  of  the  form: 


u(x,0)  =  Up  =  known 


T  (x ,  0)  =  T *  =  known 


initial  conditions  ; 


u(0,t)  = 


essential  I  u(L,t)  = 

boundary  < 
conditions  j  T(0,t)  i 


u^  =  known  or  P(o,t) 


ut  *  known  or  P(L,t) 


=  known  or  q(0,t) 


P®  =  known 


=  PL  =  known  j  natural 


T(L,t) 


Tt  =  known  or  q(L,t) 


=  q°  =  known 
=  q^  =  known 


boundary 
conditions.  (10) 


It  is  now  assumed  tho.t  a  =  a(x)  so  that  equation  (2)  reduces  to 


P  =  OA 


Therefore,  substituting  (4)  into  (5)  and  this  result  into  (11)  gives 


P  =  EA 
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The  above  result  is  now  substituted  into  (1)  to  obtain 
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which  represents  the  differential  equation  relating  displacements  and  temper¬ 
ature  to  the  applied  load  p  (x)  . 

X 

Equations  (4)  and  (7)  are  next  substituted  into  energy  balance  law  (8) 
and  this  result  is  integrated  over  the  cross-sectional  area  A  to  obtain 


where  it  has  been  assumed  that  all  field  variables  depend  on  x  and  t  only. 
Now  define 


(15) 


Careful  inspection  of  equations  (13)  and  (14)  will  indicate  that  these 
equations,  together  with  internal  state  variable  growth  laws  (6)  and  initial 
and  boundary  conditions  (9)  and  (10),  represent  a  well-posed  boundary  value  prob¬ 
lem  in  terms  of  the  2+z  dependent  variables  u,  T,  and  a  . 
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Solution  Procedure 


The  field  problem  is  to  be  solved  analytically  using  the  semi-discretized 
finite  element  technique  with  timewise  finite  differencing.  In  order  to  ac¬ 
complish  this,  differential  equations  (13)  and  (14)  must  first  be  written  in 
a  suitable  variational  form. 

Variational  Equations 

Consider  first  equation  (13),  This  governing  equation  is  integrated  against 

a  suitably  smooth  test  function  v  =  v(x)  over  the  domain  of  some  element  : 

e 

x^  <  x  <  x  . ,  : 
e  e+1 
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(16) 
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dx  =  0 


Integrating  by  parts  results  in 
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Substituting  equation  (12)  into  the  boundary  term  thus  results  in 
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Now  consider  equation  (14).  Once  again  the  governing  equation  is  inte¬ 
grated  against  a  suitably  smooth  test  function  w  =  w(x)  over  the  domain  of  the 

element  : 
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Integrating  the  heat  flux  term  by  parts  results  in 


where  equation  (15)  has  been  substituted  into  the  boundary  terms. 


Finite  Element  Spacial  Discretization 

Quadratic  displacement  and  linear  temperature  fields  are  now  chosen  within 
each  element: 


u(x,t)  =  I  u® 
i=l 


x  <  x  <  x  ..  ,  and 

e  e+1 


(21) 


T(x,t)  *  I  T* 
i=l 


x  <  x  <  x 
e  e+1 


(22) 


where  u*  -  u^(t)  and  =  T^(t)  are  the  nodal  displacements  and  temperatures, 
respectively,  and  ■  ij^(x)  and  =  <£®(x)  are  quadratic  and  linear  shape 
functions,  respectively  [17]. 

Appropriately,  v  and  w  are  endowed  with  the  properties  of  u  and  X,  respec¬ 


tively,  so  that 


(23) 


v  =  i  -  1.2,3 
w  =  4>2  i  -  1,2 


Time  dependence  in  equations  (6)  and  (20)  Is  handled  via  finite  differ¬ 
encing.  Although  higher  order  approximations  may  be  used,  Euler  forward  dif- 
ference  approximations  are  now  entered  for  the  time  rate  of  change  of  a^,  T^, 

and  ue . 
m 

3cx* 

(x,t)  =  [a®  (x,t  +  At)  -  o£  (x,t)]/At,  k  =  1,  ...,z  (24) 

dTe 

-r—  (t)  s  [T*  (t  +  At)  -  Te  (t) ] /At,  m  =  1 ,2  ,  and  (25) 

at  m  m 

due 

-r-2  (t)  s  [u*  (t  +  At)  -  u*  (t) ]/At,  m  =  1,2,3.  (26) 

at  mm 


Substitution  of  equations  (21)  through  (26)  into  equations  (18)  and  (20) 
will  result  in  (See  Appendix.) 


5x5 

where  all  nonlinearity  is  contained  in  [S],  { Fe } ,  and  {F®},  and  all  terms  are 
as  defined  in  the  appendix. 


P2  +  Pl+1  “  0  *  and  (28) 

4>2  +  -  0  .  (29) 

Boundary  conditions  are  implemented  in  the  standard  way:  1)  essential 
boundary  conditions  are  handled  by  placing  one  on  the  diagonal  of  the 

appropriate  row  and  zeros  off  diagonal  in  the  stiffness  matrix,  and  the  speci¬ 
fied  value  of  the  essential  variable  on  the  right  hand  side;  and  2)  natural 
boundary  conditions  are  implemented  directly  to  the  right  hand  side. 


Solution  of  the  Nonlinear  Algebraic  System 


Initial  conditions  are  used  for  the  first  time  step.  The  time  step  At 


is  supplied  for  each  load  increment  and  boundary  conditions  are  incremented 
directly  from  supplied  input  functions. 


The  internal  state  variable  is  handled  in  equations  (A8)  and(A22)  by 
using  equations  (24).  is  initialized  according  to  reference  18.  The  non¬ 
linear  stiffness  matrix  [S]  is  initialized  using  nodal  temperatures  and  displace¬ 
ments  from  the  previous  time  step.  The  displacements  and  temperatures  at  time 
t+At  are  then  estimated  directly  and  without  iteration  by  utilizing  equations  (27) 
for  very  small  time  steps. 
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EXAMPLE  PROBLEMS 


In  order  to  completely  define  an  example  problem  it  is  necessary 
to  specify  internal  state  variable  growth  laws  (6).  Numerous  models 
have  been  proposed  for  crystalline  metals  [18,19].  Since  it  is  not  the 
purpose  of  this  research  to  compare  these  models,  a  relatively 
established  model  proposed  by  Bodner  and  Partom  [20]  has  been  chosen. 
This  model  contains  two  internal  state  variables:  the  inelastic  strain 
(a^)  and  the  drag  stress  (a2).  The  growth  laws  for  these  variables  are 
given  by 

'  _  _2_  d  0 

a  i  '  VT  o  exP 

and 

a  2  -  m(Zi  -  a2)  odi-  A1Zi 

where  D^,  n,  m,  Z  Z^,  and  r  are  experimentally  determined  material 
constants. 

For  the  purpose  of  modeling  the  temperature  gradient  in  a  specific 
component,  a  hypothetical  problem  has  been  chosen  using  material 
properties  representati ve  of  Inconel  100  at  1005°K  (1350°F).  The 

material  and  geometric  properties  are  given  in  Table  I.  The  geometry  is 
representative  of  a  cylindrical  uniaxial  bar  which  is  2.50  inches  long 
and  0.25  inches  in  diameter.  It  has  previously  been  shown  that  Bodner 
and  Partom’ s  model  accurately  predicts  the  stress-strain  behavior  of 
INI  00  under  uniaxial  loading  conditions  for  both  monotonic  and  cyclic 
strain  controlled  loadings  [12,18]. 
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Utilizing  the  material  properties  described 


above,  th 


following  effects  have  been  studied  using  the  model  developed 
herein : 

1)  the  effect  of  variation  of  strain  rate  on  the  stress- 
strain  behavior  of  a  monotonical 1 y  extended  uniaxial  bar  which  is 
insulated  on  the  longitudinal  boundaries  (Fig.  1); 

2)  the  effect  of  strain  rate  on  the  time  dependence  of 
temperature  at  the  midpoint  of  the  bar  described  in  1)  (Figs.  2- 
*0 ; 

3)  the  spacial  temperature  variation  for  the  case  described 
above  ( Fig .  5 ) ;  and 

*0  the  effect  of  end  temperature  boundary  conditions  on  the 
temperature  at  the  center  of  a  uniaxial  bar  which  is  held  at 
fixed  temperature  at  the  end  points  and  subjected  to  cyclically 
imposed  end  displacements  (Figs.  6  and  7). 

The  slight  instability  shown  at  the  lowest  strain  rate  in  Fig. 
1  is  numerical  rather  than  an  actual  material  instability.  This 
may  be  corrected  by  using  slightly  smaller  time  steps.  However, 
this  was  not  done  herein  because  of  the  massive  computation  time 
required  with  the  current  algorithm.  Furthermore,  this  numerical 
instability  has  little  effect  on  the  predicted  temperature  field. 

The  opposite  signs  for  temperature  change  in  tension  and 
compression  shown  in  Figs.  2  through  4  is  caused  by  the  well- 
known  bulk  deformation  effect  of  thermoelastic  coupling,  which  is 
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described  by  the  second  derivative  on  displacement  in  equation 
(14).  Of  course,  after  yield,  both  tensile  and  compressive 
loadings  cause  heating  due  to  entropy  generation  caused  by  the 
inelastic  strain,  which  is  described  by  the  leading  term  in 
equation  (14). 

It  is  found  from  the  analytical  results  that  if  a  specimen 
is  mounted  in  an  experimental  apparatus  which  has  massive  grips 
simulated  by  a  fixed  temperature  boundary  condition  there  can  be 
substantial  axial  temperature  gradients  induced  in  a  time 
dependent  boundary  layer  near  the  ends  of  the  specimen.  As  shown 
in  Fig.  5,  this  boundary  layer  occurs  over  a  small  region  near 
the  end  of  the  bar  for  moderate  strain  rates  and  for  the  material 
considered  herein.  These  boundary  conditions  do  not  appear  to 
substantially  affect  the  predicted  stress-strain  behavior  (Fig. 
1),  especially  when  the  strain  measurement  is  taken  between  the 
thermal  boundary  layers  near  the  grips.  Therefore,  it  would 
appear  that  the  standard  procedure  for  obtaining  stresses  and 
strains  in  uniaxial  bars  is  not  substantially  affected  by 
mechanically  induced  axial  temperature  gradients  when  the  grips 
are  at  fixed  temperature  equivalent  to  the  initial  specimen 
temperature  and  the  bar  is  loaded  monoton i ca 1  1  y.  However,  it 
should  be  noted  that  massive  grips  which  are  mounted  outside  a 
furnace  could,  by  their  much  lower  temperature  than  the  initial 
specimen  temperature,  induce  significant  error  in  predicted 
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strains  if  the  strain  is  measured  by  dividing  relative  displacement  by  some 
gage  length. 

The  final  example  demonstrates  that  under  cyclic  loading  conditions  the 
above  conclusions  may  not  necessarily  be  true,  especially  when  the  specimen  is 
subjected  to  high-cycle  fatigue  and  at  high  strain  rates.  There  is  definitely 
a  trend  towards  an  increasing  mean  temperature  in  the  bar,  and  this  mean  tempera¬ 
ture  is  strongly  affected  by  the  thermal  boundary  conditions  as  well  as  the 
loading  rate.  Although  it  would  be  interesting  to  determine  the  mean  temperature 
rise  in  a  cyclic  fatigue  test,  the  current  algorithm  precludes  this  analysis  due  to 
the  extremely  large  computer  times  necessary  to  predict  only  a  few  cycles  of 
response  (approximately  43.8  CPU  minutes  on  an  Amdahl  470/V6  for  the  example 
demonstrated  in  Figs.  6  and  7). 

Example  3  also  demonstrates  another  interesting  phenomenon  which  may  be 
significant  in  large  space  structures.  If  the  bar  is  perfectly  insulated  the 
mean  temperature  rise  per  cycle  for  the  relatively  slow  loading  rate  shown  in 
Fig.  6  is  3.7°K,  whereas  if  the  ends  of  the  bar  are  held  at  a  fixed  temperature  of 
1005°K,  the  mean  rise  is  1.0°K  per  cycle.  Faster  loading  rates  show  less  difference 
between  the  adiabatic  result  and  the  fixed  end  temperature  result.  Since  many 
of  these  structures  are  expected  to  be  extremely  flexible  truss-like  configurations, 
a  typical  metallic  member  which  undergoes  some  yielding  (which  might  be  desirable 
in  order  to  induce  natural  damping)  might  in  fact  undergo  substantial  enough 
heating  during  vibrational  response  such  that  the  material  properties  could  be 
adversely  affected,  thus  resulting  in  a  material  related  failure  of  the  structure. 
However,  further  investigation  is  needed  on  this  last  issue  since  it  is  expected 
that  the  primary  form  of  heat  flux  off  of  space  structures  will  be  via  radiation 
on  the  longitudinal  surfaces  of  the  truss  member.  Since  the  current  analysis 


CONCLUSION 


The  current  research  has  attemped  to  demonstrate  the  effects  of 
mechanical  loading  on  one-dimensional  temperature  gradients  in  a  class 
of  viscoplastic  media.  Due  to  the  nonlinearity  and  stiffness  of  the 
field  equations,  it  was  necessary  to  utilize  a  numerical  algorithm. 
This  algorithm  has  been  shown  to  be  very  inefficient  for  solving  even 
one-dimensional  examples.  Therefore,  it  is  apparent  that  significant 
refinement  of  the  procedure  will  be  necessary  before  multi-dimensional 
analyses  can  be  performed  by  this  method.  Specifically,  it  would  be 
significant  to  determine  the  effect  of  transverse  temperature  gradients 
on  the  stress-strain  behavior  of  constitutive  specimens.  Furthermore, 
the  effects  of  thermal  boundary  conditions  on  the  longitudinal  surface 
needs  attention.  The  author  is  currently  studying  a  perturbation 
technique  for  more  efficient  solution  of  these  issues. 

The  above  points  notwithstanding,  the  current  research  demonstrates 
some  important  results.  These  are: 

1)  The  axial  temperature  gradient  in  a  viscoplastic  uniaxial  ba 
is  strongly  affected  by  the  thermal  boundary  conditions  on  the  ends. 

2)  The  end  temperature  boundary  conditions  can  cause  temperature 
gradients  which  are  substantial  enough  to  induce  spacial  variations  in 
stress  and  strains  which  invalidate  the  standard  procedure  of  using 
average  quantities,  although  when  grips  are  mounted  within  a  furnace  at 
spacially  constant  temperature,  it  appears  that  the  standard  procedure 
is  accurate. 
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3)  There  is  a  trend  toward  increasing  average  temperature  in 
cyclically  loaded  bars;  whether  or  not  this  effect  is  significant  is 
strongly  dependent  on  the  thermal  boundary  conditions  and  the  loading 
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APPENDIX 


Substitution  of  equations  (12)  and  (21)  through  (23)  into  variational  principle 
(18)  results  in 
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The  above  may  be  written  in  the  following  compact  form 
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Similarly,  substitution  of  equations  (21)  through  (23)  into  equation  (20)  re¬ 
sults  in 
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Equations  (A6)  may  be  written  in  the  following  form: 
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Time  dependence  in  equations  (6)  and  (A7)  is  handled  via  finite  differ¬ 
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Substitution  of  (A10)  through  (A12) into  finite  element  equations  (A7)  gives 
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Equations  (A14)  may  be  written  equivalently  as  follows: 
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The  above  equations  may  be  adjoined  with  equations  (A2)  to  obtain  the  following 
set  of  nonlinear  equations  for  each  element. 
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where  all  nonlinearity  is  contained  in  [S],  {F6},  and  {F6}.  Equations  (A26) 
are  identical  to  equations  (27)  in  the  main  text. 
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Abstract 


This  paper  focuses  on  the  effect  of 
thermomechanlcally  induced  heating  on  the 
response  of  a  single  member  of  a  space  truss 
structure  which  behaves  vlscoplastleally.  The 
governing  equations  are  given  for  a  typical  truss 
member,  wherein  material  inelasticity  is 
reflected  in  constitutive  equations  via  a  set  of 
internal  state  variables,  each  characterized  by  a 
history  dependent  growth  law.  The  governing 
equations  are  coupled  in  the  sense  that 
temperature  and  displacement  are  dependent  on 
each  other.  This  difficulty,  together  with  the 
fact  that  the  inelastic  constitutive  equations 
are  nonlinear  and  numerically  stiff,  requires 
that  a  computationally  complex  semldlscretized 
finite  element  spatial  technique  be  utilized  to 
obtain  a  solution.  This  procedure,  detailed 
herein,  is  utilized  to  predict  the  response  of  a 
typical  metallic  space  truss  member  under 
vibrational  or  cyclic  loading.  Particular 
interest  is  placed  on  the  temperature  rise  in 
such  a  member  due  to  hysteretic  loss  durir, 
structural  vibrations  and  in  the  presence  of 
complex  thermal  boundary  conditions 
representative  of  space  conditions.  Example 
cases  are  constructed  for  a  typical  cylindrical 
bar  of  6061-T6  aluminum  both  with  and  without 
special  coatings.  Results  indicate  that 

significant,  possibly  even  catastrophic,  heating 
can  occur  due  to  thermomechanical  coupling. 


Nomenclature 


t  -  time 

P  -  axial  Internal  resultant  force 


Px  -  axial  externally  applied  force  per  unit 


length 

x  -  axial  coordinate  dimension 
0  -  axial  stress  component 

A  -  cross-sectional  area 


T^  -  end  traction  in  units  of  force  per  unit  area 
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surface  area 

area  of  the  longitudinal  surface  of  the  bar 

■  axial  strain  component 

-  axial  displacement  component 

•  Internal  state  variable  representing  axial 
inelastic  strain 

•  Young’s  modulus  in  the  axial  coordinate 
direction 

■  coefficient  of  thermal  expansion  in  the 
axial  coordinate  direction 

■  temperature 

•  reference  temperature  at  which  no 
deformation  is  observed  at  zero 
load 

•  Internal  state  variable  representing  drag 
stress 

heat  flux  vector 
axial  component  of  heat  flux 
coefficient  of  axial  thermal  conductivity 

■  specific  heat  at  constant  elastic  strain 
mass  density 

■  internal  heat  source  per  unit  mass 

1  length  of  the  structural  element 


V 


1.  m,  Z1  ,  Zj,  ZQ,  r  -  material  constant  jsed 


in  Bodner  and  Partom’s  model 

-  flux  on  longitudinal  boundary 

-  absorbing  portion  of  perimeter  of  an  element 
normal  to  longitudinal  axis 

-  solar  radiation  flux 

-  earth  radiation  flux 
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aa  -  absorptivity 
Fg  -  earth  radiation  view  factor 
*3  -  Incident  angle  of  solar  radiation  on 
structural  component 

1E  -  incident  angle  of  earth  radiation  on 
structural  component 

o  -  Stefan-Boltzmann  constant  -  5.775  x  10 
s 

MPa  m/sec/ ( “K)** 

Tq  -  deep  space  temperature 


Introduction 


It  Is  well  known  that  In  viscoplastic  metals 
a  certain  amount  of  mechanically  Induced 
hysteretlc  mechanic  energy  loss  Is  converted  to 
heat,  thus  resulting  in  a  temperature  rise  In  the 
medium.  In  recent  research* ’ J  a  model  has  been 
developed  ror  predicting  this  effect  by  utilizing 
thermodymanlc  constraints  together  with 
constitutive  equations  of  Internal  state  variable 
type*.  Furthermore,  It  has  been  shown  that  in  a 
perfectly  Insulated  uniaxial  bar  ,  as  well  as  In 
a  uniaxial  bar  with  Insulated  longitudinal 
surface  and  fixed  end  temperature'*,  significant 
temperature  rise  can  occur  In  the  component 
during  cyclic  loading. 

The  purpose  of  the  current  research  Is  to 
simulate  the  response  of  a  typical  metallic  space 
truss  structural  element  (see  Fig.  1 )  In  the 
postylelded  state  and  to  determine  If  significant 
heating  occurs  when  this  component  Is  subjected 
to  cyclic  mechanical  loading.  This  problem  Is  of 
Interest  because  a  certain  amount  of  material 
Inelasticity  is  desirable  In  order  to  produce 
passive  structural  damping.  The  factors  of 
interest  In  this  simulation  are  the  effects  of 
thermal  boundary  conditions  and  loading  rate  on 
the  thermal  response.  In  particular,  It  is  of 
interest  to  determine  If  radiative  boundary 
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Fig.  1.  Typical  Space  Truss  Structural  Element. 
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conditions  on  the  longitudinal  surface  of  the 
truss  component  are  significant  enough  to  carry 
off  all  heat  generated  due  to  hysteretlc  loss. 

The  paper  first  reviews  the  governing  field 
equations,  then  briefly  discusses  the  procedure 
whereby  a  numerical  algorithm  Is  constructed  for 
modeling  the  problem.  This  Is  followed  by  a 

detailed  discussion  of  the  Implementation  of 
thermal  boundary  conditions.  Finally,  example 
results  are  obtained  for  representative  space 
structural  components. 
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Governing  Field  Equations 

The  governing  field  equations  were  presented 
in  reference  5  for  quasi-static  conditions.  For 
problems  Involving  Inertial  effects,  the 
governing  equations  are  as  follows: 

a)  equilibrium^’ 

H.  •  p,(x) 

3x  * 

(1) 

where  the  axial  resultant  P  Is  defined  by 

P  -  /odA 
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b)  strain-displacement  relation 
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c)  thermomechanical  constitution, 
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where  a,  and  a-  are  the  internal  state  variables 
( ISV)  representing  inelastic  strain  and  drag 
stress,  respectively,  In  the  constitutive  model 
developed  by  Bodner  and  Partom  .  Several  other 
constitutive  models  have  been  developed  for 
viscoplastic  metals,  and  these  are  reviewed  in 
references  1  and  8.  Finally, 

d)  conservation  of  energy^, 

[(Ee-Ea.+EaT0)3al  +  Ea2T3T]-EaT3e  -PCST-iSl+er  *° 
1  R  at  3 t  3t  3t  3x 
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(9) 

Conservation  of  mass  is  satisfied  trivially 
(under  the  assumption  of  small  motions  in  a 
closed  system),  and  the  second  law  of 
thermodynamics  has  been  prevloi^^y  shown  to  be 
satisfied  by  the  above  equations  ’  . 

The  governing  equations  are  adjolnei^  with 
appropriate  initial  and  boundary  conditions  such 
that  a  well-posed  boundary  value  problem  is 
constructed  In  terms  of  the  following  dependent 
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variables  which  are  sought  as  functions  of  x  and 
o.c.u.q.T.P.m, .  and  a,.  Due  to  ISV  growth 


t:  a.c.u.q.T.P.a, ,  and  a,.  Due  to  ISV  growtn 
laws  (6)  and  (7)  (as  well  as  radiative  boundary 
conditions),  the  problem  is  nonlinear. 


Solution  Procedure 


As  described  in  detail  in  reference  5  for 
the  quasistatic  problem,  the  solution  is  obtained 
using  the  semi-discretized  finite  element 
technique,  wherein  finite  elements  are 
constructed  spatially,  and  finite  differencing  is 
iised  in  time.  The  result  is  a  time  marching 
algorithm  which  is  reviewed  here  briefly. 

First,  equations  (4)  and  (5)  are  substituted 
into  (2)  and  this  result  is  substituted  into  (1) 
to  give  the  following  equllbrium  equation: 


{ EA  f  3u  - 

ara(T-TR)1  ) 

l  La* 

Ji 

-  Px(x) 


(10) 


Next,  equation  (4)  is  substituted  into  equation 
(9)  to  obtain  the  coupled  energy  balance  law. 


-EaT32u 


3t3x 


V.q  +pr  -  0 


(11) 


-PC  3T 
v  — — 

3t 

The  result  is  a  set  of  two  coupled  partial 


differential  equations  in  terms  of  axial 
displacement  u«u(x,t)  and  temperature  T-T(x.t). 


Variational  Principles  and  Finite  Element 


Discretization 


Selecting  a  suitably  smooth  test  function 
v«v(x )  over  the  domain  of  some  element  8  : 
x_<x<x_,,  one  may  construct  the  following 
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var 
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lational  principle  from  equation  (12)  : 


•  f  6+1  EA  3v  r3u-a.-a(T-T  )1  dx 

Jx  1  A. 


-v(x_,)P(x 


e+r  e+1 


.Xpl.! 

)+v(x  )P(x vp  dx 


(12) 


where  the  boundary  terms  result  from  the  standard 
integration  by  parts. 

The  variational  principle  for  heat  equation 
(11)  is  constructed  by  first  integrating  this 


equation  against  a  test  function  w-w(x)  on  8g  to 


obtain 


W  (Ee-Ea^>EaT^)^al-t-Ea2T  3T  J 


-EaT  3  u 


3x3t 


-pC  3T  +  7-q  +  pr 
v  3t 


dV 


(13) 


Integrating  the  flux  term  by  parts,  assuming  that 
nonaxlal  components  of  flux  are  negligible,  and 
substituting  equation  (8)  will  thus  result  in 
x 


I  w  I A  fE3u-Ea +EaT„\  3al  +  Ea2T  3T 

i  III"  )~  >'J 

e  2  } 

-AEaT  3  u  -  ApC  3T  -  kA  3w  £T  ( dx 
3t3x  vTt"  3x  3x  I 

/•xe+ 1 

“'u(xe+l)Aq(xe+l)  +  w(xe)Acl(%)-  I  c  “  9cd* 


dx 


(14) 


Variational  equations  (12)  and  (14)  are  now 
discretized  by  assuming  the  following 
displacement  and  temperature  fields  in  a  typical 
element  (superscripted  e): 


1  e  e 

u(x.t)-  tmluL  (t)^t(x)  xe<  x<xg+1 


T(x,t)-£  T*  (t)*®(x)  xg<x<xe+1 


(15) 

(16) 


where  ,.e 

ui 

temperatures 


and  T® 


are 


quadratic 


respectively, 
and  linear 


nodal  displacements 
and  \(i®  and  $e 


and 

are 


1  i 

shape  functions, 
respectively'.  Furthermore,  v  and  w  are  endowed 
with  the  properties  of  u  and  T.  Note  that  a 
higher  order  element  must  be  used  for 
displacement  than  temperature  due  to  the  fact 
that  temperature  produces  strain  rather  than 
displacement. 

Timewise  discretization  is  Implemented  via 
the  following  backward 
equations: 


dT  (t)s[T„(t)-T  (t-At)]Mt 
q  mm 


dt 


du*  (t)*[„!(t)-u®(t-6t)]Mt 

m  ®  ® 

dt 

The  above  equations  require 


finite 

difference 

u»-l,2 

(17) 

m-1.2.3 

(18) 

small  time  steps  in 
order  to  guarantee  numerical  accuracy.  However, 
they  are  uncondltonally  stable  which  is  necessary 
because  ISV  gr^gth  laws  (6)  and  (7)  are 
numerically  stiff  . 

Substitution  of  equations  (15)  through  (18) 
into  the  governing  field  equations  in  variational 
form  will  result  in  the  following  algebraic 
equations: 


3x3  i  3x21 
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(19) 


5x1 


5x1 


[K°], 


[se]. 


[Ke], 


where 

described  in  reference  5,  and 


Cse], 


and 


(Fe) 


F?s?e 


f*e+l 

l  i  < 

J  * 


dx 


(20) 


where  is  as  defined  in  reference  5. 


The  last 

term  in' the  above  equation  accounts  for  thermal 
flux  boundary  conditions  on  the  longitudinal 
surface  of  an  element. 

After  global  assembly  and  imposition  of 
boundary  conditions  equations  (19)  can  be  solved 
in  a  time  marching  scheme  in  order  to  obtain  the 
nodal  displacements  and  temperatures  as  functions 
of  time. 

Global  assembly  of  the  element  equations  is 
accomplished  ^n  the  standard  way  using  the 


Boolean  matrix 


Imposition  of  Boundary  Conditions 


For  a  typical  space  truss  structural 
element,  the  boundary  conditions  are  assumed  to 
be  of  the  following  type: 


u(0,t)  •  u“  -  known 


u(L,t)  *  u^  «  known 


(21  ) 


T(0 ,  t ) 


known 


T(L,t)  *  T  -  known 

t 


:  V  V  V.V.  V  \\V.V  V.V.  .\\\ 


qc“'as[q3  C0S  X3+FE(1'aE)q3C09Ac+FEqEcosAe] 


where  the  first  term  in  the  above  equation  is  the 
solar  radiation  flux  absorbed  by  the  body,  the 
second  term  is  the  solar  radiation  flux  reflected 
by  the  earth  and  absorbed  by  the  body,  the  third 
term  is  the  earth  radiation  flux  absorbed  by  the 
body,  and  the  last  term  is  the  flux  radiated  by 
the  member  to  space. 

The  above  boundary  condltons  may  be 
implemented  to  the  <Usoretlzed  global  equations 
in  the  standard  way  .  Although  equation  (22) 
technically  Includes  the  unknown  temperature 
field,  the  component  temperature  is  treated  as  a 
known  quantity  in  this  term  for  each  time  3tep. 
This  approximation  is  acceptable  due  to  the  fact 
that  the  numerical  stiffness  of  constitutive 
equations  (6)  and  (7)  requires  extremely  small 
time  steps  in  order  to  obtain  an  accurate 
solution. 

EXAMPLE  PROBLEMS 

A  typical  structural  element  has  been 
modeled  with  properties  shown  In  TABLE  1  .  The 
material  properties  were  obtained  experimentally 
In  the  Mechanics  and  Materials  Center  at  Texas 
A&M  University  for  A1  5086  at  room  temperature, 
which  1s  similar  to  A1  6061-T6. 

Sample  cases  were  constructed  for  various 
cyclic  loading  rates  for  two  different  sets  of 
thermal  boundary  conditions,  as  described  in 
TABLE  2.  Both  cases  are  considered  to  be  "worst 
cases"  in  that  the  component  is  In  a  maximum 
radiation  flux  condition  at  the  maximum 
equilibrium  temperature  during  one  orbital  cycle. 
The  two  cases  differ  In  the  emmlssivtty  and 
absorptivity  values  for  the  component  due  to 
differences  in  surface  treatment  of  the 
component.  For  case  I,  the  component  Is 
anodized,  and  for  case  II,  the  component  is 
painted  with  high  emmlssivlty  ITTRE-S1 3GL0  white 
paint  . 

We  now  consider  two  elements  In  a  large 
space  structure  (see  Fig.  1).  Both  elements  are 
constructed  of  the  3ame  material  and  are 
geometrically  identical.  However,  element  one  is 
painted  with  the  high  emmlssivlty  paint  described 
above  and  is  in  full  view  of  both  earth  and  sun, 
whereas  element  two  Is  anodized  and  13  In  view  of 
earth  only.  For  thi3  case,  as  described  in  Table 
2,  the  components  have  Identical  equilibrium 
temperatures  Tp  -  295°K(obtained  by  setting  q  *0 
in  equation  (22*)  )• 

In  both  cases  the  structural  members  have 
been  subjected  to  50  cycles  of  loading  at  three 
different  frequencies:  1  Hz,  5  Hz,  and  25  Hz. 
These  frequencies  have  been  selected  as 
representative  of  resonant  frequencies  in  a 
representative  space  structure.  For  example,  a 
typical  structure  analyzed  in  reference  14  has 
resonant  frequencies  of  4.1  Hz  and  J.4  Hz  in  the 
first  two  modes.  Because  the  resonant  frequency 
of  the  first  mode  in  the  structural  element 
itself  is  240  Hz,  inertial  effects  may  be 
neglected  in  these  examples. 

Results  for  the  cases  described  above  are 
shown  in  Figs.  2  through  8.  In  Figs.  2  through  4 
the  cyclic  stress-strain  curve  is  shown  at  the 


900  J/kg/°K  (0.215  Btu/lb/ °F) 
23.8x10'6in./ln./°K  ( 1 3- 2x1 06  in./ln./°F) 
1.27  x  10""  MPa  m2/sec/°K 
(73. 4  Btu/ft/h/°F) 

71.0  x  IQ3  MPa  (10.3  x  106  psl) 

6.45  x  10'"  m2  (1 .00  in2) 

295  °x  (72°F) 

3.66  m  (12.0  FT) 


10  x  10  m/m 
1.685  x  10"7  sec  _1 


0.1770  MPa'1  (1.2205  Ksl'1 ) 


620.1  MPa  (89.93  Ksl) 


2.66  Mg/m3  (0.096  lb/in.3) 
0.0508  m  (0.8333  Ft.) 
387.8  MPa  (56.25  Ksl) 


TABLE  1 .  Material  and  Geometric  Properties  for  a 
Typical  Truss  Structural  Element  (from  reference 
11). 

CASE  I  CASE  II 
a  0.20  (degraded)  0.3218  (degraded) 


1 . 39  MPa  m/sec 


0.20  MPa  m/sec  0.20MPa  m/sec 


( 4, 080  km  altitude) 


296. 2°K  ( 73 ■ 6°F) 


296. 2°K 


CASE  I  -  Surface  painted  with  S13GL0  white 


CASE  II  -  Chromic  anodized  surface 


T'toLE  2.  Thermal  Properties  for  Example  Cases  I 
and  II  (from  references  12  and  13). 


location  x-L/2  for  CASE  I  and  at  all  three 
loading  rates.  It  is  found  that  in  all  cases  the 
specimen  reaches  cyclic  saturation  after 
approximately  five  cycles.  Thereafter,  the 
hysteretlc  energy  loss  per  cycle  becomes  a 
constant  value. 

In  Figs.  5  through  7  the  temperature  rise  is 
plotted  for  both  cases  at  all  three  loading 
rates.  As  expected,  the  amount  of  temperature 
rise  Increases  dramatically  with  loading  rate. 
For  example,  after  50  cycles  the  total 
temperature  rise  at  x-L/2  is  17.5 °K£ 1  Hz), 

62. 5°K(5  Hz),  and  119. 7°K  (25  Hz)  for  case  I. 
Furthermore,  it  is  apparent  that  while  neither' 
surface  treatment  can  be  regarded  as  resulting  in 
negligible  heating;  at  the  higher  loading  rates 
the  anodized  surface  treatment  produces 
temperature  rises  which  are  significantly  higher 
than  those  where  the  surface  is  painted  with 
ITTRE-S1 3GL0  paint.  Finally,  it  is  believed  by 
these  researchers  that  the  nonlinear  nature  of 
the  average  temperature  rise  per  cycle  suggests 
that  the  temperature  rise  asymptotically 
approaches  some  upper  bound,  although  this  belief 
cannot  be  corroborated  at  this  time  due  to  the 
large  computer  times  required  in  the  current 
algorithm. 

Fig.  8  shows  that  the  spatial  temperature 
variation  at  5  Hz  is  approximately  spatially 
homogeneous.  Apparently,  a  very  thin  boundary 
layer  forms  near  the  end  of  the  component,  and 
this  boundary  layer  has  little  effect  on  the 
temperature  at  x-L/2.  In  fact,  subsequent 
investigations  by  the  authors  have  shown  that,  at 
least  for  the  geometry  and  physical  conditions 
considered  herein,  identical  results  may  be 
obtained  more  efficiently  by  neglecting  spatial 
variations  in  displacement  and  temperature. 
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Fig.  3.  Cyclic  Stress-Strain  Curve  at  x-L/2  for 
Case  I  Coating  Loaded  at  5  Hz. 
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Fig.  H.  Cyclic  Stress-Strain  Curve  at  x-L/2  for 
Fig.  2.  Cyclic  Stress-Strain  Curve  at  x-L/2  for  Case  I  Coating  Loaded  at  25  Hz. 

Case  I  Coating  Loaded  at  1  Hz. 
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Conclusion 


The  currant  research  has  attempted  to 
predict  the  response  of  a  typical  space 
structural  element  which  Is  viscoplastic  and  is 
subjected  to  various  cyclic  loading  conditions  in 
the  presence  of  radiation  boundary  conditions. 
Several  general  conclusions  can  be  made  as  a 
result  of  this  reasearch: 

1)  significant  temperature  rises  may  occur 
due  to  hysteretic  loss,  although  the  precise 
amount  depends  on  loading  rate  and  surface 
treatment; 

2)  the  special  paint  ITTRE-S13GL0  appears  to 
produce  significantly  lower  temperature  rises 
than  anodized  surface  treatment; 

3)  the  temperature  rise  appears  to  be 
approaching  an  upper  bound  which  Is  loading  rate 
and  surface  treatment  dependent;  and 

*1)  the  thermal  boundary  layer  which  forms 
near  the  end  of  the  member  appears  to  have  little 
effect  on  the  far-fleld  temperature  rise. 

These  conclusions  Indicate  that  future 
research  on  this  subject  should  perhaps 
concentrate  on  spatial  variations  In  the  radial 
direction  rather  than  the  axial  direction.  More 
Importantly,  these  results  Indicate  that  an 
inelastic  structural  component  may  undergo 
temperature  rises  during  structural  vibrations 
which  are  so  substantial  that  the  material 
properties  of  the  component  may  be  further 
degraded,  thus  leading  to  rallure  of  the 
component  and  perhaps  even  failure  of  the  entire 
structure. 
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Abstract 


In  this  paper  the  effect  of  degradation  of 
material  properties  on  structural  frequencies  and 
mode  shapes  of  Large  Space  Structures  (LSS)  Is 
Investigated.  The  difficulty  and  cost  of 
maintenance  of  LSS  make  It  a  necessity  to  design 
these  structures  to  operate  with  a  certain  amount 
of  load-induced  damage.  This  damage  Is  commonly 
observed  in  fibrous  composite  media. 


Sensitivity  studies  conducted  on 
representative  space  truss  structures  indicate  that 
degradation  of  material  properties  may  have  a 
significant  effect  on  the  structural  mode  shapes 
and  frequencies.  For  even  small  amounts  of 
reduction  In  stiffness  (10J),  frequencies  and  nodal 
locations  may  change  significantly.  It  is  clear 
that  these  effects  must  be  taken  Into  consideration 
when  designing  control  systems  for  Large  Space 
Structures. 


Introduction 


Due  to  economic  constraints,  It  Is  projected 
that  advanced  high  strength-to-welght  ratio 
aerospace  materials  will  be  utilized  in  future 
generation  space  structures.  Such  materials 
include  polymer  and  metal  matrix  fibrous 
composites,  which  are  known  to  .undergo  a  certain 
amount  of  load  induced  damage.  ’  These  materials 
are  also  expected  to  undergo  a  certain  amount  of 
environmentally  Induced  damage  or  degradation,  thus 
resulting  in  significant  stiffness  losses. 


Experimental  research  on  advanced  composite 
materials  indicates  that  the  material  may  undergo 
up  to  15  percent  loss  In  stiffness  due  to 
thermomechanical  fatigue,  which  causes  a  variety  of 
damage  modes  in  the  structure.  Additional  loss  of 
stiffness  may  be  attributed  to  elevated  temperature 
and  chemical  changes  due  to  solar  radiation  and 
other  environmental  effects.  Thl3  reduction  In 
stiffness  affect3  the  dynamic  response  which  In 
turn  1 3  critical  In  the  development  of  control 
systems  Tor  LSS.  In  this  paper,  sensitivity 
studies  will  be  presented  which  Investigate  the 
effect  of  stiffness  loss  on  structural  frequencies 
and  mode  shapes. 


The  advent  of  the  space  shuttle  has  made 
possible  the  development  of  LSS.  Control  systems 
for  stabilizing  and  maneuvering  these  very  large 
space  structures,  especially  those  for  precise 
pointing,  will  require  extension  of  current 
technology. 
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Although  large  size  by  Itself  does  not  arouse 
concern,  structural  flexibility  resulting  from 
minimizing  the  structural  weight  In  non- 
gravitatlonal  fields  may  present  problems. 
Extremely  large  structural  flexibility  may  result 
in  large  amplitudes  and  low  frequencies  (.01  to  10 
Hz)  which  may  create  new  complications  for  control 
designers. 


J 


As  an  example  of  the  precision  required  ,  a 
typical  radlometry  application  may  utilize  a  200 
meter  antenna  with  an  effective  beam  width  of  0.01 
degrees  and  have  requirements  limiting  the 
vibratory  beam  shift  to  le33  than  0.005  degrees  and 
dynamic  surface  distortions  to  less  than  1mm. 
Maneuvering  or  maintaining  the  altitude  of  such  a 
satellite  leads  to  flexible  body  motion  which  mu3t 
be  well  predicted  and  controlled. 


The  importance  of  interaction  between  control 
systems  and  vibratory  response  has  causejJ 
considerable  research  in  LSS  control  systems. 

The  current  practice  of  guaranteeing  a  large 
separation  between  modal  frequencies  and  the 
bandwidth  of  control  will  not  be  adequate  in  future 
applications.  The  combination  of  large  size  and 
payload-weight  restrictions  will  drive  structural 
frequencies  down  and  the  need  for  more  accurate 
pointing  will  drive  the  control  system  bandwidth 
up.  When  sufficient  frequency  separation  becomes 
impossible,  there  exists  a  need  Tor  adaptive 
control  systems.  This  leads  to  further  research  in 
the  design  of  structural  control  systems  actuator/ 
sensor  placement,  and  distributed  sensing  and 
actuation  as  opposed  to  co-located  sensors  and 
actuators. 


Techniques  for  achieving  modal  control  of  LSS 
will  require  a  more  accurate  knowledge  of  modal 
characteristics.  Optimum  sensor  and  actuator 
placement  will  be  greatly  Influenced  by  modal 
effects  which  must  be  known  to  a  greater  degree  of 
precision. 


Problem  Summary 


In  order  to  investigate  the  possible  effects 
of  material  degradation  on  the  dynamic  response  of 
LSS,  a  representative  space  truss  structure  nas 
been  selected  In  the  shape  of  a  long  boom  as  shown 
in  Fig.  1.  Using  several  loading  histories,  stress 
distributions  have  been  obtained  for  each  truss 
member .  The  resulting  stress  distributions  can  be 
used  in  a  material  damage  model  to  define  material 
degradation  and  resultant  stiffness  reductions. 
Using  the  reduced  stiffness  properties,  modal 
analyses  have  been  conducted  on  the  structure  to 
show  the  effect  of  material  degradation  on  natural 
frequencies,  mode  shapes  and  nodes.  Details  of  the 
finite  element  model,  material  degradation  model, 
and  numerical  results  are  presented  below. 
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Fig.  1  Space  Truss  Structure 


Model  Description 


Material  Degradation  Model 


The  process  of  ultimate  failure  of  composite 
materials  is  preceded  by  a  sequence  of 
microstructural  and  macro3tructural  events  which 
are  termed  as  damage.  These  events  may  be  due  to 
transverse  cracking,  delamlqatJUsn,  fiber  breaking 
and  fiber-matrix  debonding.  The  mechanical 

response  of  the  structure  is  affected  by  this 
damage.  Global  material  properties  like  stiffness 
and  residual  strength  may  be  substantially  alteped 
during  the  life  of  the  structural  components. 
Some  of  the  analytical  studies  tj<g-  modeling  damage 
Include  ^  shear  lag  concept,  fracture  based 

concepts,  .  and  internal  state  variable 


concepts , 1  ^  and  internal  state  variable 

theories.  Although  Important  progress  has 


been  made,  current  understanding  of  damage  is  not 
complete. 


Damage  in  polymeric  composites  is  modelled  in 
this  paper  a3  a  load  history-dependent  reduction  in 
stiffness  in  each  structural  element.  The  internal 
state  variable  theory  (ISV)  is  used  for  modeling 


mechanical 


relationship  13  of  the  form, 


°ij  '  C’ijkl  (tkl  *  ekl  5 


In  this  case,  the  ISV  are  assumed  to  be 
second  order  tensor  valued  and  to  enter  only 


through  the  modulus  tensor. 


effective  modulus  tensor  given  by 


P  mnkl IJ  a  mn  ’  P"1 • • • • «r 


are  a  set  of  r  internal  state 


variables  which  are  given  by  the  following  set  of 
ISV  growth  laws. 


nPmn(Ekl’T,akiq) 


At  low  homologous  temperatures  these 
materials  are  assumed  to  be  rate  insensitive  so 
that  the  above  model  will  result  in  quasi-elastic 
(rate  independent)  equations  in  which  inelasticity 


is  reflected  only  through  the  slowly  degrading 
modulus  tensor.  Experimental  evidence0' 


Indicates  that  the  time  scale  for  degradation  of 
C-  .  is  very  long  compared  to  the  frequencies  and 
mode  shapes  of  representative  structures.  It  Is 
therefore  sufficient  for  many  space  structural 
applications  to  treat  equations  0)  in  the  degraded 
state  only. 


The  stress-strain  relationship  for  the  truss 
elements  is  a  one-dimensional  approximation  of 
equations  (2)  given  by 


°xx  -  E’(£xx  -  e  xx> 


where  <>xx  and  exx  are  the  uniaxial  stress  and 


strain,  e  is  the  thermal  strain,  and  E’  is  the 
axial  stlfWiess  of  the  truss  element  given  by 


E  (1  -  a) 


where  E  is  the  undegraded  axial  stiffness  and  a  is 
a  scalar  valued  parameter  representing  the 
integrated  effect  of  all  damage  modes  such  as 
matrix  cracking,  interlaminar  fracture,  fiber 
breakage,  and  fiber-matrix  debonding. 


Experimental  research  on  composite  matsrlals 
Indicates  a  power  law  degradation  of  axial 

(  ffnaaa  a«  a  f  1  An  rtf  cfroqq  hi  efrtMu  *  ' 


stiffness  as  a  function  of  stress  history.  * 
Hence  the  damage  ISV  growth  law  is  assumed  to  be  of 
the  form 


k1  ( o/omax) 


where  k1  and  n  are  material  parameters,  o  is  the 
maximum  stress  in  the  structure,  and  o  is  the  axial 


stress  in  each  truss  element.  For  constant  stress 
amplitude,  equation  (6)  may  be  Integrated  in  time 
to  give  the  following  approximation 


aCV  -  k',  [o(t,)/omaxr 


where  k’  and  n1  are  material  parameters  which  may 
be  time  dependent. 


A  power  law  form  of  damage  is  used  herein  for 
simplicity  and  for  an  Initial  attempt  at  modeling 
the  structural  response  with  damage.  In  reality 
the  damage  laws  will  be  more  complex  1  and  are 
currently  being  developed  for  future  work. 


Finite  Element  Model 


Figure  1  lllus-rates  the  geometry  of  the 
representative  space  truss  used  to  simulate  an 
antenna  boom.  This  structure  is  sixty  feet  long 
with  10  bays,  six  feet  long  by  three  feet  wide. 
The  finite  element  model  has  ’24  space  truss 
elements  and  44  nodes.  In  the  Initial  undegraded 
configuration,  the  material  properties  are  the  same 
for  all  members  with  the  following  values: 


2 


Material  type:  Graphite  epoxy  (Hexel) 

Young's  modulus  E  -  21.5x10°  gal 
Cross  sectional  area  ■  1.0  In 
Density  *  0.065  lb/in^ 

Coefficient  of  thermal  expansion  -  2x10  0  ln/in/° F 
Reference  temperature  -  89.6°F 

Each  truss  member  is  idealized  with  a 
standard  six  degree  of  freedom  truss  element  of 
constant  cross  section.  Because  the  structure  is 
idealized  as  linear  with  slowly  varying  material 
properties,  conventional  linear  finite  element 
methodology  may  be  used  t^.wnlte  global  equations 
of  equilibrium  of  the  form'  ' 

CM]{q}  *  [K]{q)  -  (Q)  (8) 

where  [M]  Is  the  mass  matrix,  [K]  Is  the  stiffness 
matrix,  { q }  is  the  nodal  displacement  vector,  and 
(Q)  Is  the  nodal  force  vector.  The  stiffness 
matrix  [K]  is  dependent  on  the  spatially  variable 
damage  state  a  which  varies  from  element  to 
element.  Standard  eigenvalue  extraction  may  be 
performed;  in  this  case,  subspace  Iteration  was 
used  to  obtain  the  first  five  frequencies  and  mode 
shapes. 


The  spatial  distribution  of  degradation  and 
stiffness  reduction  of  LSS  will  be  complex  and 
dependent  on  loading  and  environmental  history. 
For  the  present  Investigation,  wherein  material 
degradation  is  assumed  to  be  a  function  of  stress 
history.  It  wa3  necessary  to  make  some  assumptions 
about  the  corresponding  stress  history  and  spatial 
distribution  of  stresses  within  the  L.SS. 

Two  approaches  were  used  to  obtain  candidate 
stress  histories/di3tributions  for  predicting  the 
stirfness  degradation.  In  one  case,  the  stress 
distribution  was  obtained  for  an  assumed  thermal 
load  history/dlstrlbutlon.  Secondly,  a  modal 
approach  was  used  wherein  it  was  assumed  that 
primary  degradation  occurred  in  the  first  two 
bending  mode3  of  the  structure.  After  computing 
the  mode  shapes  for  the  first  two  undegraded 
bending  modes,  the  nodal  displacements  were  used  to 
compute  a  corresponding  stress  distribution. 

In  each  case,  the  degradation  model  given  by 
equation  (7)  was  then  U3ed  to  obtain  degraded 
properties  for  each  trus3  member  assuming  that  the 
element  stressed  the  most  was  degraded  a  specified 
percentage.  The  resultant  structure  with  degraded 
properties  has  spatially  variable  stiffness  that 
varies  from  element  to  element.  Mode  shapes  and 
frequencies  were  then  computed  with  varying  maximum 
percent  ig°3  of  degraded  properties. 


Discussion  of  Results 

Natural  frequency  and  mode  shape  responses 
have  been  obtained  for  several  stress-induced 
degradation  te3t  cases  as  described  above  for  the 
"•epresentat  1  ve  space  truss  structure  shown  in  Fig. 
1.  Thi3  particular  truss  structure  geometry, 
representing  a  segment  of  a  boom,  is  similar  to 
ones  being  used  for  other  PACOSS  related  work. 
Assuming  the  boom  is  fixed  on  one  end  (at  x»0),  the 
five  lowest  frequencies  (for  the  virgin  structure) 
are  equal  to  3 - *•  Hz,  9.5  Hz,  9.6  Hz,  19.2  Hz,  and 
20.3  Hz.  The  first  mode  is  a  combined  torslon- 
inplane  shear  mode,  the  next  two  modes  are  bending 


modes  about  the  z  and  y  axes,  respectively,  and  the 
fourth  mode  Is  a  pure  torsion  mode. 

The  first  case  considers  the  boom  structure 
shown  in  Fig.  1  (which  is  assumed  to  be  fixed  on 
one  end)  with  a  thermal  gradient  over  the  cross- 
section.  It  is  likely  that  one  surface  of  the 
space  structure  will  become  significantly  hotter 
than  the  other  surface  due  to  solar  heating, 
attitude  of  the  structural  elements  and  shadowing 
effects.  To  Investigate  the  effect  of  this  thermal 
gradient  through  the  depth  of  the  truss,  the 
stresses  In  each  element  were  calculated  by 
specifying  a  temperature  of  12?°F  for  the  members 
on  the  top  surface,  80.6°F  for  th“  ...embers  on  the 
bottom  surface  and  100°F  for  the  diagonal  members 
connecting  the  top  and  bottom  surface.  With  this 
thermally-induced  stress  distribution,  the  axial 
stiffness  of  each  element  was  degraded  by  using 
equation  (7).  The  maximum  level  of  degradation 
(loss  of  stiffness)  was  set  to  a  prescribed 
percentage  for  the  element  with  the  highest  stress 
and  remaining  elements  were  degraded  according  to 
their  stress  level  by  using  equation  (7).  The 
value  of  n’  In  equation  (7)  was  assumed  to  be  0.75. 

In  Fig.  2  the  first  three  natural  frequencies 
are  plotted  for  different  levels  of  damage.  The 
effect  of  damage  on  the  natural  frequencies  is 
clear.  Increasing  the  level  of  damage  reduces  the 
stiffness  of  the  space  truss  and  this  In  turn 
drives  the  natural  frequencies  down  significantly 
even  for  modest  damage  states.  For  a  maximum  loss 
of  25X  In  axial  stiffness  (for  the  highest  stressed 
members),  the  first  three  natural  frequencies  are 
reduced  by  about  8%.  Since  mode  shapes  are 
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Fig.  2  Effect  of  Damage  on  Natural  Frequencies 
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important  for  designing  the  control  systems  of  the 
large  space  structures,  it  is  desirable  that  they 
be  constant  with  time,  Although  it  was  found  that 
there  was  no  appreciable  change  in  the  first  mode 
shape  between  the  undegraded  and  degraded  cases, 
higher  modes  were  altered  due  to  material 
degradation.  Figure  3  is  a  plot  of  the  z 
displacement  for  the  second  mode  shape  along  the 
length  of  the  space  truss  (z-0,  y«0).  Significant 
changes  in  the  mode  shape  and  node  locations  as  a 
function  of  percent  degradation  are  observed.  The 
sign  of  the  modal  displacement  is  reversed  near  the 
free  edge  for  the  degraded  and  undegraded  cases  and 
the  location  of  the  node  (zero  displacement) 
changes  appreciably.  Figure  1  is  :  similar  plot  of 
the  y  displacement  along  the  length  of  the  space 
truss  for  the  third  mode. 
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The  value  of  n'  in  equation  (7)  was  varied 
from  0.25  to  1.0  to  study  its  effect  on  the  mode 
shapes.  It  was  found  that  the  trend  in  mode  shape 
changes  was  similar  for  different  n' .  Figure  5 
illustrates  this  point.  Here  the  z  displacement 
for  the  second  mode  is  plotted  along  the  length  of 
the  space  truss  for  different  values  of  n'  (maximum 
reduction  in  axial  stiffness  was  20t).  The  plot 
indicates  that  increasing  n'  (l.e.,  decreasing  the 
nonlinearity  of  the  degradation  model)  tends  to 
Increase  the  changes  In  the  modal  displacements. 
Such  nonlinearity  becomes  increasingly  Important 
when  stresses  vary  spatially  over  the  structure, 
i.e.,  some  members  are  highly  stressed  compared  to 
others. 


ig.  3  Effect  of  Degradation  on  Second  Mode 
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Fig.  4  Effect  of  Degradation  on  Third  Mode 


Effect  of  Material  Degradation  Exponent 
on  Second  Mode  Shape 


The  next  two  sample  oases  consider  the 
situation  where  we  assume  that  primary  degradation 
occurs  in  the  first  two  bending  modes.  For 
simplicity,  it  Is  assumed  that  damage  occurring  in 
one  mode  does  not  affect  the  damage  in  any  others, 
l.e.,  no  damage  Induced  coupling  of  modes.  In 
reality,  this  may  not  be  the  case  and  will  be 
considered  In  future  research. 

In  the  first  case,  we  consider  the  case  where 
degradation  has  occurred  In  the  first  bending  mode, 
l.e.,  degradation  Is  based  on  stresses  calculated 
from  the  modal  displacements  corresponding  to  the 
second  mode  shape.  Figure  6  shows  the  resulting 
first  three  natural  frequencies  for  different 
levels  of  damage.  For  a  maximum  reduction  In 
stiffness  of  25*  the  first  three  natural 
frequencies  decrease  by  8.6*,  9.2*  and  7.6*, 
respectively.  There  Is  little  change  In  the  first 
mode  shape  for  the  degraded  and  undegraded  cases. 
Figure  7  is  a  plot  of  the  z  displacement  for  the 
second  mode  shape  along  the  length  of  the  space 
truss  and  shows  that  the  modal  displacements  change 
quite  drastically  for  the  degraded  structure.  The 
displacement  at  the  free  edge  Is  nearly  30  times 
the  magnitude  of  the  undegraded  case  for  a  maximum 
damage  of  25*  (the  sign  of  the  displacement  Is  also 
reversed)  and  the  location  of  nodes  also  change 
considerably.  Figure  8  Indicates  similar  changes 
In  the  y  displacement  for  the  third  mode  shape. 
The  fourth  mode  (torsional)  Is  relatively 
unaffected  by  the  degradation  of  material  stiffness 
properties.  This  Is  as  expected  because  the 
present  analysis  assumed  that  primary  degradation 
occurred  In  the  first  two  bending  modes.  Different 
results  would  be  expected  lr  significant  stiffness 
reduction  occurred  in  the  primary  torsion  mode. 


F.g.  5  Effect  of  Damage  on  Natural  Frequencies 


Fig.  7  Effect  of  Degradation  on  Second  Mode 
Assuming  Second  Mode  Damage  State 


Fig.  8  Effect  or  Degradation  on  Third  Mode 
Assuming  Second  Mode  Damage  State 
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Results  have  also  been  obtained  for  the  case 
where  damage  Is  assumed  to  occur  In  the  third  mode 
(second  bending  mode).  As  in  the  previous  examples 
there  Is  no  appreciable  change  in  the  first  mode 
shape  between  the  undegraded  and  degraded  cases . 
The  z  displacement  corresponding  to  the  second  mode 
shape  13  plotted  in  Fig.  9  for  different  levels  of 
damage.  The  displacement  at  the  free  edge  Is  very 
large  In  the  damaged  states  as  compared  to  the 
undegraded  state.  Figure  10  Illustrates  similar 
results  for  the  third  mode  shape.  These  results 
show  that  the  mode  shapes  and  node  points  may 
change  significantly  for  even  small  damage  amounts. 


Conclusions 

This  study  has  attempted  to  Investigate  the 
possible  effects  of  material  damage  and  stiffness 
reduction  on  the  modal  response  of  LSS.  Large 
space  structures  constructed  of  fibrous  composites 
will  experience  some  stiffness  reductions  produced 
by  load-induced  and  environmentally-induced  damage 
of  the  material.  To  what  extent  this  will  occur  Is 
uncertain  at  this  point  but  even  small  damage 
amounts  appear  to  be  significant. 

The  present  work  has  shown  that  load-induced 
degradation  of  material  properties  may  have  a 
significant  effect  on  the  structural  frequencies 
and  mode  shapes.  For  the  representative  boom 
structure  considered  here,  even  small  amounts  of 
material  stiffness  degradation  (10J)  produce 
frequency  and  node  shifts  which  appear  to  be 
significant.  It  is  not  Inconceivable  that  mode 
shapes,  node  locations,  and  frequency  distributions 
will  change  over  the  plant  design  life  in  such  a 
way  that  the  structure  response  Is  very  much 
different  from  the  virgin  structure.  Such  changes 
in  plane  response  would  require  "robust"  control  Of 
3  nature  which  may  not  be  possible  with  present 
technology.  Consequently,  it  is  Important  that 
these  effects  be  taken  into  consideration  when 
designing  the  control  systems  for  large  space 
structures . 

Although  preliminary,  this  study  suggests  the 
need  f or  a  more  accurate  knowledge  of  the  physical 
nature  of  material  degradation  in  fibrous 
composites,  it3  influence  on  structure  stiffness, 
and  how  material  degradation  will  affect  the  long¬ 
term  modal  characteristics  for  large  space 
structures . 
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ABSTRACT 


A  general  constitutive  theory  is  formulated  for  a  short 
fiber  metal  matrix  composite.  The  constitutive  theory  is 
based  on  continuum  mechanics  with  constraints  provided  by 
fracture  mechanics  and  thermodynamics.  The  basic  premise 
of  the  constitutive  theory  is  that  load-induced 
microstructural  damage  associated  with  the  inclusion  of  the 
fibers  in  the  metal  matrix  results  in  a  loss  of  material 
st i f f  ness . 

The  concept  of  the  damage  dependent  constitutive  theory 
was  evaluated  by  an  experimental  investigation  of  the 
microstructural  damage  in  an  aluminum  6061-T6  matrix  with 
silicon  carbide  ( SiC )  whiskers.  The  primary  objective  of 
the  investigation  was  the  identification  of  microstructural 
damage  such  as  voids  or  cracks  that  were  associated  with 
the  presence  of  the  SiC  whiskers. 

The  results  of  the  experimental  investigation  include 
measured  reductions  in  the  elastic  modulus  of  tensile 
specimens  which  were  loaded  beyond  yield  and  unloaded. 

Also,  there  was  an  associated  change  in  the  extent  of  the 
microstructural  damage.  These  results  lead  to  the 
conclusion  that  the  effect  of  microstructural  damage  must 
be  included  in  the  constitutive  relationship  in  order  to 
fully  describe  the  behavior  of  a  metal  matrix  composite. 

KEY  WORDS 

Composite  material,  metal  matrix,  damage,  fatigue, 


ites,  the  continuous  fibers  in  the  metal  matrix  composites 
can  be  oriented  in  different  directions  to  tailor  direc¬ 
tional  strength  and  stiffness  to  satisfy  particular  design 
requirements . 

To  effectively  utilize  the  capabilities  of  metal 
matrix  composites,  a  model  that  simulates  the  behavior  of 
the  material  must  be  developed.  Since  composites  develop 
load  induced  microstructural  fracture  during  their  operat¬ 
ional  life,  this  damage  must  be  incorporated  into  the  model 
in  order  to  construct  a  valid  constitutive  model.  Both 
theoretical  and  experimental  efforts  are  required  to 
develop  the  model. 

This  paper  presents  a  model  that  will  predict  the 
thermomechanical  constitutive  behavior  of  randomly  oriented 
short  fiber  quas i - i so tr op ic  metal  matrix  composites  under 
both  monotonic  and  cyclic  fatigue  loading  conditions.  The 
model  incorporates  various  effects  such  as  damage, 
inelastic  strain,  and  viscoelasticity  in  the  constitutive 
equations  with  constraints  obtained  from  continuum 
mechanics  and  thermodynamics  with  internal  state  variables. 

The  concept  of  the  damage  dependent  constitutive  model 
requires  experimental  verification.  Also,  the  full 
development  of  the  model  requires  an  experimentally 
generated  data  base.  Therefore,  this  paper  also  presents 
the  results  of  an  experimental  program  aimed  at  addressing 
the  above  two  issues.  The  material  system  selected  for 
study  is  a  6061-T6  aluminum  matrix  with  chopped  silicon 
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carbide  whiskers.  The  objective  of  the  experimental 
program  is  to  determine  if  there  is  a  load  induced 
reduction  in  the  elastic  modulus  of  tensile  test  specimens 
along  with  an  associated  observable  physical  phenomenon 
such  as  a  change  in  the  state  of  microstructural  damage. 

PREVIOUS  RESEARCH 

A  vast  amount  of  research  has  been  reported  on  com¬ 
posite  materials.  Most  of  the  work  involves  experimental 
studies  on  the  behavioral  response  of  materials.  To  a 
lesser  extent,  general  constitutive  models  have  been  de¬ 
veloped.  General  information  pertaining  to  metal  matrix 
composites  is  presented  by  Vinson  [1]  and  Renton  [2],  These 
authors  discuss  fabrication  techniques  for  both  continuous 
fiber  and  particulate  reinforced  metal  matrix  composites  as 
well  as  selected  applications  and  properties  of  the 
material.  Fabrication  techniques  include  powder 
metallurgy,  liquid  metal  infiltration,  diffusion  bonding, 
electroforming,  rolling,  extrusion,  pneumatic  impaction  and 
plasma  spray.  There  have  been  a  number  of  investigations 
C3“19]  that  have  studied  the  physical  behavior  of  metal 
matrix  composites.  Divencha,  Fishman,  and  Karmarkar  [7] 
have  studied  aluminum  with  short  whisker  (SiC) 
reinforcement.  The  authors  found  that  the  Al/SiC  composite 
properties  varied  significantly  with  the  amount  of  SiC 
whiskers  used.  Maclean  and  Misra  [8]  have  examined  the 
applications  of  Al/SiC  composites  for  aerospace  use.  They 
found  that  although  the  strength  and  stiffness  of  the 
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material  were  higher  than  that  of  conventional  aluminum, 
fracture  toughness  decreased. 

In  order  to  observe  damage  in  metal  matrix  composites, 
various  experimental  techniques  have  been  investigated. 
Pipes,  Ballintyn,  Scott,  and  Carlyle  [10]  have  observed  the 
acoustic  emission  response  of  various  laminates  of  B/Al 
composites.  Both  mechanical  and  acoustic  emission 
responses  of  the  composites  are  presented.  Various 
non-destructive  examining  techniques  such  as  ultrasonic 
C-scan,  radiography,  acoustic  emissions,  stiffness  change, 
edge  replication,  and  eddy  current  were  investigated  by 
Ulman  and  Henneke  [11]  to  evaluate  large  scale  imperfect¬ 
ions  and  damage.  X-ray  radiography  and  ultrasonic  C-scan 
were  found  to  be  excellent  methods  to  observe  both  initial 
imperfections  and  damage  due  to  loading. 

Numerous  investigative  groups  [12-19]  have  reported 
the  observation  of  damage  in  metal  matrix  composites. 
Several  different  material  systems  have  been  studied, 
including  B/Al ,  BSiC/Ti,  and  Al/SiC.  It  has  been  found 
that  complex  damage  states  develop  in  both  continuous  riber 
laminated  composites  and  chopped  fiber  composites. 

A  number  of  theoretical  papers  have  also  been  published 
on  constitutive  equations  for  metal  matrix  composites 
[20-26],  However,  to  date  these  efforts  have  not  included 
a  damage  parameter.  Rate  dependent  viscoplastic 
constitutive  models  have  been  developed  for  neat  metals 
[27~39].  These  models  have  not  included  a  damage  parameter 


applicable  to  composites.  Theoretical  models  have  been 
developed  for  predicting  constitutive  properties  of  cracked 
bodies  [40-50].  However,  to  date  these  models  have  not  been 
applied  to  metal  matrix  composites.  A  review  of  the  current 
literature  indicates  that  no  constitutive  model  has  yet 
been  developed  which  is  applicable  to  metal  matrix 
composites  that  undergo  large  scale  matrix  yielding  and 
load  induced  damage. 

MODEL  DEVELOPMENT 

This  section  presents  a  short  review  of  a  constitutive 
theory  for  composite  media  with  damage  which  is  described 
in  detail  in  references  49  through  52.  The  model  is 
presented  herein  because  it  provides  the  motivation  for  the 
experimental  research  program  undertaken  as  part  of  this 
research. 

The  constitutive  framework  is  based  on  a  continuum 
mechanics  approach  with  constraints  on  the  relations 
provided  by  thermodynamics  and  fracture  mechanics.  The 
general  model  is  applicable  to  materials  with  damage  (such 
as  voids,  cracks,  etc.)  and  includes  inelastic  effects  such 
as  plasticity  [52].  The  model  is  constructed  within  the 
framework  of  continuum  mechanics  and  thermodynamics.  The 
governing  conservation  laws  are  integrated  over  a  small 
local  volume  element  which  is  assumed  to  have  a 
statistically  homogeneous  damage  state,  as  shown  in  Fig.  1. 
The  Helmholtz  free  energy  can  be  expressed  as 
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Helmholtz  free  energy  due  to  the  elastic-plastic  response 
in  the  absence  of  damage,  and  u£  i3  the  energy  due  to 
damage.  It  is  therefore  hypothesized  that 


hEP  -  hEP(£ij ,ejj ,AT) 


(2) 


where  e  ^  is  the  total  strain  tensor,  e^  is  the  inelastic 
strain  tensor,  and  AT  is  the  temperature  difference  from 
the  reference  temperature.  Furthermore, 


UL  "  UL( eij ’ eij ’ AT,aij ) 


(3) 


where  a,,  is  the  internal  state  variable  representing 
damage,  which  is  defined  by 


(V 


where  S  is  the  surface  area  of  cracks  in  the  local 
volume  VT  ,  and  u.  and  n  are  the  crack  opening 
displacements  and  normals,  as  described  in  Fig.  2. 
Constraints  imposed  by  the  second  law  of  thermodynamics 
give  the  following  result  [51]: 
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Therefore,  expanding  equations  (2)  and  (3)  in  Taylor  series 
expansions  in  terms  of  their  arguments,  substituting  into 
equation  (5)  and  truncating  higher  order  terms  results  in 


1  i j  "  °ij  +  Cijkl(ekl  "  ekl  ‘  °kl  ‘  ekl) 


(6) 
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where  Ojj  is  the  residual  stress  tensor,  is  the 

inelastic  strain  tensor,  is  the  thermal  strain  tensor, 

and  is  the  linear  elastic  modulus  tensor.  For  the 

uniaxial  case  in  which  there  is  negligible  temperature 

change,  the  following  form  results: 


o  =  oR  +  E^(e  '  -  a)  (7) 

where  EL  is  the  initial  loading  elastic  modulus.  Now 
define  the  initial  unloading  modulus  Ey  such  that  (See  Fig. 
3)  : 
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It  is  assumed  that  at  relatively  low  homologous 
temperatures  the  inelastic  strain  remains  constant  when 


unloading,  so  that: 


Assuming  linearly  elastic  unloading  of  the  matrix,  the 


change  in  damage  is  proportional  to  the  change  in  strain: 
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constant  (upon  unloading)  -  8 


(10) 


Therefore  for  unloading 
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The  model  described  above  provides  the  motivation  for 
the  experimental  research.  It  is  hypothesized  that  the  two 


parameters  that  are  defined  in  this  model  (a^  and  8)  can 


be  determined  by  experimental  methods.  By  determining  the 
change  between  the  initial  loading  and  unloading  moduli  of 
the  composite  in  a  uniaxial  mechanical  test,  8  (defined  in 


equation  11)  can  be  found.  It  is  also  observed  that 


(defined  in  equation  4)  can  be  determined  by  evaluating  the 
amount  of  surface  area  in  the  composite.  It  is  therefore 
desirable  to  determine  if  a  cause-and-ef f ect  relationship 


exists  between  the  microstruc tural  damage  (a,  .  )  and  the 


stiffness  loss  ( g) 


EXPERIMENTAL  PROGRAM 


As  stated  before,  the  primary  objective  of  the 
experimental  effort  is  to  provide  support  and  documentation 
for  the  constitutive  model  for  metal  matrix  composites 
developed  in  the  previous  section.  This  initial 
experimental  study  will  provide  an  indication  of  the  type 
of  damage  that  may  be  present  in  the  "virgin"  material  and 
an  indication  of  any  additional  damage  that  is  created  from 
mechanical  loading  of  the  specimen.  The  term  damage 
includes  cracks,  voids,  and  debonding  of  particles  and 
matrix  material.  All  of  these  modes  of  damage  cause  new 
surface  area  to  be  created  in  the  material.  One  of  the 
goals  of  this  effort  is  to  identify  and  quantify  the  amount 
of  damage  in  the  composite.  To  accomplish  this  objective, 
it  is  necessary  to  develop  an  experimental  evaluation 
technique  and  mechanical  test  program.  The  procedure  must 
be  capable  of  inducing  and  then  examining  microstructural 
damage  in  the  composite. 

Material  and  Specimen  Fabrication 

The  material  used  in  this  study  was  obtained  from  ABCO 
Metals  Silag  Operation  in  Greer,  S.  C.  The  composition  of 
the  material  is  6061  aluminum  with  a  twenty  percent  volume 
fraction  of  F-9  grade  silicon  carbide  whiskers.  The  SiC 
whiskers  average  two  microns  in  diameter  and  twenty  microns 
long.  The  constituents  were  cast  into  a  billet  form  by  a 
powder  metallurgy  process  and  then  rolled  into  a  plate.  The 
composite  has  a  T-6  temper.  Prior  to  specimen  machining, 
the  plate  was  examined  by  C-scan  to  observe  whether  or  not 
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were  computed  and  compared  at  several  load  levels. 
Mlcrostructural  Damage  Evaluation 

Once  the  specimens  were  mechanically  loaded,  it  was 
necessary  to  examine  the  microstructure  for  damage.  It  was 
desirable  to  use  a  nondestructive  technique  in  the 
evaluation  process  since  the  most  convenient  way  to  observe 
the  initiation  and  growth  of  damage  is  to  study  the 
behavior  of  a  particular  region  of  one  specimen.  However, 
the  photographic  resolution  of  NDE  techniques  such  as 
ultrasonics  or  x-ray  radiography  was  not  adequate  to 
distinguish  the  mlcrostructural  damage  of  interest,  which 
could  be  as  small  as  one  tenth  of  a  micron.  Therefore,  it 
was  decided  that  scanning  electron  microscopy  ( SEM )  would 
be  the  best  method  of  mlcrostructural  evaluation. 

Resolution  of  the  photomicrographs  obtained  from  SEM  is 
better  than  one  tenth  of  a  micron.  The  disadvantage  of 
using  SEM  is  that  the  technique  is  destructive  by  nature. 
This  implies  that  growth  of  damage  in  one  particular 
specimen  under  various  conditions  cannot  be  measured. 

The  preparation  of  the  internal  surfaces  of  the  tensile 
specimens  followed  a  procedure  that  was  established  after 
evaluating  the  application  of  several  methods  [53, 5^]  to 
the  A1/S1C  specimens.  The  specimens  were  sectioned  into 
pieces  about  one-half  by  one-half  inch  (12.5  mm  X  12.5  mm) 
for  mounting  in  the  scanning  electron  microscope.  Since 
there  was  uncertainty  about  the  degree  of  anisotropy  of  the 
fiber  orientation  in  the  plane  of  the  rolled  plate,  two 
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sections  were  cut  from  each  coupon  so  that  two  orthogonal 
views  of  each  coupon  could  be  examined.  The  sectioned 
pieces  were  then  mounted  in  a  conductive  mounting  material 
(Konductomet  I)  manufactured  by  Buehler.  After  mounting, 
the  specimens  were  polished  with  diamond  paste,  chemically 
etched,  and  cleaned  in  an  ultrasonic  cleaner.  It  should  be 
noted  that  specimen  surface  preparation  requires  great  care 
because  of  the  vast  difference  between  the  hardness  of  SiC 
and  aluminum.  Unless  care  is  taken  and  an  appropriate  grade 
and  hardness  of  abrasive  compound  is  used,  an  uneven 
terrain  may  be  created  by  the  removal  of  the  softer 
aluminum  matrix  leaving  exposed  SiC  whiskers. 


The  final  step  in  specimen  preparation  was  to  vacuum 
deposit  a  thin  (about  300  angstroms)  coating  of  gold  to  the 
viewing  surfaces.  This  is  necessary  to  achieve  a  good 
image  on  the  SEM.  When  the  voltage  is  applied  to  the 
specimens  by  the  SEM  ,  free  electrons  are  released  from  the 
specimens.  It  is  these  electrons  that  are  received  by  a 
sensor  which  forms  the  specimen  image.  Light  molecular 
weight  compounds  have  less  free  electrons,  thus  the  image 
formed  is  not  as  sharp  as  an  image  of  a  compound  with  a 
large  molecular  weight.  Since  gold  is  an  element  with  a 
large  molecular  weight,  a  thin  coating  on  the  specimens 
provides  an  electron  source  without  losing  surface  detail. 

The  specimens  were  examined  in  a  Jeol  JSM-25  II 
scanning  electron  microscope.  After  choosing  the  optimum 
acceleration  voltage  for  the  type  of  specimens,  a  series  of 


SEM  photographs  were  taken  of  the  surfaces  of  orthogonal 
edges  of  the  specimens  under  both  virgin,  loaded,  and 
failed  conditions.  These  photographs  included 
magnifications  ranging  from  2.000X  to  10.000X.  Along  with 
the  photographs  of  these  polished  surfaces,  additional 
photographs  of  fracture  surfaces  were  taken  to  observe  the 
effect  of  the  whisker  re  inf orcement . 

The  photomicrographs  were  qualitatively  and 
quantitatively  evaluated  for  microstructural  damage 
associated  with  the  SiC  particles.  A  scheme  for 
determining  the  actual  surface  area  of  the  damage  was  based 
on  standard  quantitative  me tal lograph ic  techniques  [55]. 

The  evaluation  method  consisted  of  covering  the  photograph 
with  a  transparent  grid  of  twenty  by  twenty  divisions  per 
inch.  A  count  was  made  of  the  intersections  of  the  grid 
that  covered  damage.  The  percent  of  damage  surface  area 
was  then  obtained  by  dividing  the  intersections  covering 
the  damage  by  the  total  number  of  intersections  and 
multiplying  by  one  hundred. 

RESULTS 

As  noted  before,  the  results  presented  herein  were 
obtained  for  the  primary  purpose  of  evaluating  the  concept  of  the 
damage  dependent  constitutive  relationships.  No  attempt  was  made 
to  construct  a  complete  data  base  to  support  the  full  development 
of  the  constitutive  model.  A  comparison  of  the  typical 
stress-s tra i n  response  of  the  A1  6Q61-T6,  Al-SiC  with  0° 
orientation  and  Al-SiC  with  90°  orientation  is  given  in  Fig.  5. 
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As  can  be  seen,  the  effect  of  adding  the  SIC  particles  is  to 
increase  the  elastic  modulus  and  to  increase  the  yield  and 
ultimate  strengths.  These  typical  stress-strain  curves  were  very 
useful  for  determining  the  intermediate  strain  levels  used  to 
evaluate  the  unloading  moduli. 

Several  tensile  specimens  were  monotonically  loaded 
past  yield,  followed  by  continued  loading  with  periodic 
unloading  to  obtain  unloading  moduli  at  different  maximum 
load  levels.  This  test  scheme  allowed  for  the 
determination  of  the  initial  loading  modulus  and  several 
unloading  moduli  at  different  load  levels.  Typical  material 
mechanical  responses  are  shown  in  Figs.  6  and  7,  for  the  0° 
and  90°  orientations,  respectively.  Values  of  initial 
loading  moduli  are  compared  to  numerous  unloading  moduli 
for  several  specimens  in  Figs.  8  and  9  for  the  0° 
orientation  and  90°  orientation,  respectively.  The  moduli 
were  calculated  using  a  least  squares  curve  fitting 
program.  It  is  obvious  that  the  unloading  moduli  decrease 
with  increasing  strain  for  both  specimen  orientations. 

Since  the  value  of  the  modulus  was  very  dependent  on  the 
location  of  the  data  points  on  the  unloading  curve,  a  consistent 
method  of  data  point  selection  was  adhered  to  during  all  modulus 
calculations.  Three  different  sets  of  data  points  were  selected 
on  each  unloading  curve  for  each  test.  The  unloading  modulus  was 
then  calculated  from  each  set  of  data  points.  The  three  sets  of 
data  points  used  were: 

a.  the  first  six  data  points  following  load 


& 


reversal , 

b.  the  middle  ten  data  points  on  the  unloading 
curve,  and 

c.  all  the  data  points  on  the  unloading  curve. 
Three  values  of  the  initial  loading  modulus  were  also 
calculated  using  the  same  technique.  The  mean  and  standard 
deviation  of  these  data  point  sets  were  then  calculated  for 
the  Initial  loading  curve  and  each  subsequent  unloading 
curve.  By  comparing  sets  of  unloading  moduli  calculated  in 
this  manner,  consistent  trends  in  the  data  were  identified. 

The  90°  specimens  exhibited  a  clear  reduction  in  unloading 
modulus  with  increasing  applied  strain.  On  the  other  hand, 
the  unloading  moduli  of  the  0°  specimens  were  less 
sensitive  to  increasing  strain. 

Photomicrographs  obtained  from  the  SEM  were  very  useful  in 
obtaining  information  about  the  microstructure  of  the  composite. 
Fig.  10  shows  a  composite  view  of  three  mutually  orthogonal 
sections,  where  the  coordinate  axes  are  the  same  as  in  Fig.  4. 

The  orientation  and  actual  shape  of  the  whiskers  should  be 
observed.  The  question  of  whether  the  SiC  whiskers  are  randomly 
oriented  in  the  plane  that  the  Al/SiC  billet  was  rolled  was 
important  in  view  of  the  differences  in  the  0°  and  90°  specimen 
data.  It  is  clear  from  the  SEM  photographs  that  the  whiskers  are 
oriented  somewhat  in  the  principal  rolling  direction  (x),  which 
coincides  with  the  loading  direction  of  the  0°  specimens. 

It  is  necessary  to  look  at  the  microstructure  of  both 
virgin  and  failed  specimens  to  determine  if  there  are  load 
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induced  changes  in  the  state  of  microstructural  damage.  The 
damage  state  was  determined  from  five  photom icr ographs  of 
randomly  selected  edge  and  side  views  of  each  loading  state 
for  a  single  coupon.  Figs.  11  and  12  are  representative 
photomicrographs  from  virgin  and  loaded  specimens  for  the 
0°  and  90°  orientations,  respectively.  Upon  close 
examination  of  the  photomicrographs,  it  was  observed  that 
there  were  numerous  voids  and  crack-like  features  of  the 
microstructural  surfaces  of  both  the  virgin  and  loaded 
specimens.  Therefore,  the  fractographic  examinations 
attempted  to  identify  only  changes  in  the  surface  features 
between  the  virgin  and  loaded  specimens.  The  most  extensive 
change  (see  Fig.  12)  in  the  m icrostructure  was  due  to 
cracks  in  the  aluminum  matrix  extending  from  the  SiC 
particles.  There  were  no  other  extensively  observed 
load-induced  microstructural  damage  modes.  Table  1  lists 
the  amount  of  damage  that  was  measured  in  each  group  of 
photomicrographs  for  specimens  loaded  in  both  the  0°  and 
90°  orientations..  It  is  apparent  from  Figs.  11  and  12  as 
well  as  Table  1  that  there  is  little  change  in  the 
microstructure  of  damage  for  specimens  loaded  in  the  0° 
direction,  while  specimens  loaded  in  the  90°  direction 
exhibit  significant  load-induced  microstructural  damage. 


The  most  dramatic  photomicrographs  of  the 


microatructural  response  of  the  Al/SiC  material  were 
fracture  surface  photos.  One  of  the  most  important 
considerations  in  the  determination  of  crack  onset  in 
composite  materials  is  the  bond  between  the  reinforcing 
particles  and  the  matrix  material.  Fig.  13  is  a  low 
magnification  (1 50X)photograph  of  the  fracture  surface. 

When  the  fracture  surface  is  viewed  under  a  higher 
magnification  (3.000X,  see  Fig.  14)  two  major  observations 
can  be  made.  First,  the  many  cusps  in  the  matrix  show  that 
on  the  microatructural  level,  the  material  is  very  ductile. 
The  second  observation  is  that  the  SiC  whiskers  form  a  very 
strong  bond  with  the  aluminum.  This  is  due  to  the 
observation  that  the  whiskers  remain  in  the  matrix, 
although  they  appear  to  be  bonded  only  on  a  small  portion 
of  their  surface  area  (see  arrow  Fig.  25).  The  result  of 
this  strong  bond  is  that  when  stress  concentrations  exist 
at  the  interfaces  of  the  aluminum  matrix  and  reinforcing 
whiskers,  the  matrix  may  yield  and  inelast ically  deform 
before  Interfacial  cracks  initiate. 

As  stated  earlier,  specimens  with  the  0°  orientation 
(3ee  Fig.  8)  Indicated  relatively  little  load-induced 
reduction  in  modulus.  This  tends  to  suggest  that  there  is 
no  mechanism  present  to  produce  a  reduction  in  the 
stiffness  of  the  test  specimen.  This  is  supported  by  a 
comparison  of  the  photomicrographs  of  a  virgin  specimen  and 
a  loaded  specimen  (see  Fig.  11),  which  shows  no  apparent 
difference  in  microstructural  damage  state.  Furthermore,  a 
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specimen  subjected  to  1000  cycles  with  the  maximum  strain 


exceeding  yield  also  did  not  exhibit  an  increase  in  the 
microstructural  damage.  On  the  other  hand,  the  specimens 
with  the  90°  orientation  (see  Fig.  9)  exhibited  a 
significant  reduction  in  modulus  with  increasing  strain 
level.  This  trend  does  suggest  the  presence  of  a  mechanism 
that  produces  a  reduction  in  the  stiffness  of  the  test 
specimen.  This  contention  is  supported  by  the  comparison 
of  photomicrographs  of  a  virgin  specimen  and  a  loaded 
specimen  (see  Fig.  12),  which  reveals  an  obvious  and 
considerable  increase  in  the  microstructural  damage 
associated  with  the  SIC  particles.  Finally,  a  fatigue 
loaded  90°  specimen  also  exhibited  both  a  substantial 
increase  in  microstructural  damage  along  with  a  significant 
reduction  in  the  elastic  modulus.  A  similar  observation 
has  been  made  by  Johnson  in  continuous  fiber  metal  matrix 
composites  [15,16,45]. 

CONCLUSION 

A  constitutive  model  of  the  behavior  of  short  fiber 
reinforced  metal  matrix  composites  has  been  formulated.  The 
constitutive  model  was  developed  from  the  continuum  mechanics 
viewpoint  with  constraints  imposed  by  thermodynamic 
considerations.  The  model  utilizes  an  internal  state  variable  to 
characterize  the  development  of  load  induced  microstructural 
damage  associated  with  the  inclusion  of  the  fiber  particles  in 
the  metal  matrix.  The  internal  state  variable  is  defined  by  the 
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kinematic  behavior  of  the  microcracks  and  enters  the 


stress-strain  relationships  as  a  strain-like  quantity. 

It  was  shown  in  the  development  of  the  constitutive 
model  that  the  damage  parameter,  cxt j  >  can  be  used  to 
predict  stiffness  losses.  Therefore,  it  was  the  objective 
of  this  research  to  measure  stiffness  loss  and  the 
associated  microstructural  damage  as  a  function  of  strain 
level  in  order  to  qualitatively  assess  the  applicability  of 
the  model  to  the  Al-SiC  metal  matrix  composite.  This 
experimental  objective  was  carried  out  by  determining  the 
initial  loading  and  subsequent  unloading  moduli  of  tensile 
specimens  oriented  parallel  (0°)  and  perpendicular  (90°)  to 
the  principle  rolling  direction  of  a  plate  fabricated  from 
6061-T6  aluminum  with  silicon  carbide  whiskers.  In 
addition,  scanning  electron  microscopy  was  utilized  to 
characterize  and  quantify  load-induced  changes  in  the 
microstructural  damage  associated  with  the  silicon  carbide 
particles . 

The  results  showed  that  the  Al-SiC  plate  was 
anisotropic  with  approximately  1 5-20*  difference  in  the 
moduli  of  the  specimens  oriented  in  the  0°  and  90° 
directions.  Also,  the  SEM  photomicrographs  indicated  that 
the  SiC  whiskers  were  oriented  more  or  less  parallel  to  the 
principle  rolling  direction. 

There  was  very  little  difference  between  the  initial 
loading  and  subsequent  unloading  moduli  of  the  0° 
specimens.  Also,  there  was  no  apparent  load-induced  change 
In  the  state  of  microstructural  damage  in  these  specimens. 


On  the  other  hand,  there  was  a  significant  reduction  in  the 
moduli  of  specimens  oriented  in  the  90°  direction. 
Furthermore,  the  photomicrographs  revealed  very  obvious  and 
significant  load-induced  changes  in  the  state  of 
microstructural  damage  in  the  90°  specimens. 

The  results  illustrate  a  clear  cause  and  effect  between 
the  increase  in  load-induced  microstructural  damage  and  a 
decrease  in  the  elastic  modulus  of  the  Al-SiC  metal  matrix 
composite.  It  is  concluded  from  these  results  that  the 
constitutive  behavior  of  a  short-fiber  reinforced  metal 
matrix  composite  can  only  be  modelled  by  an  appropriate 
treatment  of  the  microstructural  damage  associated  with  the 
fiber  particles.  The  constitutive  model  proposed  herein 
accounts  for  this  microstructural  damage  through  a  load 
history  dependent  internal  state  variable.  Although  the 
results  must  be  considered  qualitative  at  this  time,  it 
appears  that  the  proposed  model  may  be  acceptable  for 
predicting  stiffness  loss  as  a  function  of  microstructural 
damage . 
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TABLE  1 . '  Percent  of  Damage  (Crack)  Surface  Areas 


SPECIMEN 

ORIENTATION 

SPECIMEN 

LOAD 

HISTORY 

PERCENT 

CRACK 

SURFACE  AREA* 

STANDARD 

DEVIATION 

0°  (END  VIEW) 

VIRGIN 

11.17 

2.29 

0°  (END  VIEW) 

MONOTONIC 

FAILED 

9. 48 

1  .55 

0°  (END  VIEW) 

1000  CYCLE 
FATIGUE. 

9.31 

1  .16 

0°  (SIDE  VIEW) 

VIRGIN 

9.53 

1.55 

0°  (SIDE  VIEW) 

MONOTONIC 

FAILED 

8.62 

1  .33 

0°  (SIDE  VIEW) 

1000  CYCLE 
'  FATIGUE 

8.37 

1.50 

90°  (END  VIEW) 

VIRGIN 

1  .20 

0.90 

90°  (END  VIEW) 

MONOTONIC 

FAILED 

14.7 

2.30 

90°  (END  VIEW) 

1000  CYCLE 
FATIGUE 

31  .2 

3.10 

*Both  large  voids  and  cracks  were  counted  as  damage  in  the  specimens 
with  the  0°  orientation. 
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Abstract 


Using  internal  state  variables,  the  time  dependent  behavior  of  a  thermoviscoplastic 
uniaxial  rod  is  reduced  to  the  solution  of  a  nonlinear  initial  value  problem.  Nonlinear 
black-body  radiation  boundary  conditions  of  Stefan  type  are  considered.  In  the  quasi¬ 
static  case,  one  can  solve  for  all  the  field  quantities  explicitly  in  terms  of  the  stress. 
Using  this  representation,  we  derive  upper  and  lower  bounds  on  the  temperature  vari¬ 
ation  within  the  rod.  Full  coupling  between  mechanical  and  thermal  effects  are 
allowed.  These  bounds  agree  quite  closely  with  numerical  solutions  of  the  initial 
boundary  value  problem.  Boundary  layers  form  near  the  ends  of  the  rods.  Series 
representations  for  the  compliance  and  temperature  are  derived.  Applications  to  some 
problems  involving  space  structures  are  considered.  Finally,  a  result  concerning  isoth¬ 
ermal  nonlinear  wave  dispersion  in  such  materials  is  obtained. 

INTRODUCTION 

Large  space  structures  require  significant  passive  damping  in  order  to  sus¬ 
tain  structural  integrity  during  structural  vibrations  in  a  microgravity  field. 
One  passive  damping  mechanism  which  has  been  proposed  is  material  inelasti¬ 
city.  However,  in  this  method  of  damping,  a  substantial  portion  of  the  strain 
energy  in  the  structure  is  converted  into  heat  via  hysteretic  loss.  The  purpose 
of  this  paper  is  to  derive  bounds  on  the  asymptotic  temperature  increase  in 
terms  of  bounds  on  the  cyclic  stress,  for  a  simple  model  representing  one  ele¬ 
ment  in  a  space  structure.  In  particular,  lower  bounds  are  derived  which  indi¬ 
cate  substantial  local  heating  may  take  place  if  the  member  is  operated  near  its 
limiting  strength  .  These  temperature  increases  may  in  fact  lead  to  substantial 
material  degradation  over  time. 

In  Section  7,  an  example  of  a  uniaxial  rod  coated  with  high  emissivity 
low  absorbtion  white  paint  in  an  outer  space  environment  is  considered. 
Numerical  results  are  presented  and  compared  with  the  asymptotic  bounds 
derived  in  this  paper. 
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Section  1.  Derivation  of  the  equations  of  Equilibrium.  Constitutive 
Laws. 

The  equations  we  use  to  describe  the  quasi-static  response  of  a  thermo¬ 
viscoplastic  uniaxial  bar  have  been  formulated  in  [  1  ] .  For  brevity,  we  sum¬ 
marize  the  main  points  in  this  section.  The  principal  equations  are  conserva¬ 
tion  of  momentum 


Oijj  -  0  (1.1) 

where  o, y  is  the  stress  tensor,  and  the  conservation  of  energy,  as  expressed  in 
the  modified  Fourier  heat  conduction  law  [  2  ] 

—  *kl  +  A jklteij&uTf  )  ~  (1*2) 

DijklfoijT  fid)  —  pCyT  “  CL 

DiJkl  is  the  elastic  modulus  tensor,  %  is  the  infinitesimal  elastic  strain 
tensor,  is  the  inelastic  strain  tensor,  3,;  is  the  thermal  expansion  tensor,  TR 
is  the  reference  temperature  at  which  zero  deformation  produces  no  stress,  p  is 
the  mass  density,  and  cv  is  the  specific  heat  at  constant  volume.  In  addition, 

qj  -  -kTj  (1.3) 

where  k  is  the  coefficient  of  thermal  conductivity.  Dots  indicate 

differentiation  with  respect  to  time. 

Furthermore,  we  have  the  stress-strain  relationship 

°ij  ”  A/«(««  —  *kl  ~  *L)  (1*4) 

where  ^  -  a«(T  —  T R )  is  called  the  thermal  strain  tensor. 

The  internal  state  variables  efj  obey  the  following  evolution  equations 

“  g(*U'T,a£)  (1.5) 

where 

A5  -  <*&)  (1.6) 
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are  an  additional  set  of  internal  state  variable  growth  laws,  reflecting  local 
averages  of  such  microphysical  phenomena  as  dislocation  density  and  disloca¬ 
tion  arrangement. 

We  remark  that  the  constitutive  equations  (1.4)  are  obtained  from  ther¬ 
modynamic  constraints  imposed  on  a  broad  class  of  materials  with  internal 
state  variables.  The  heat  conduction  equation  (1.2)  is  the  result  of  thermo¬ 
dynamic  constraints  imposed  on  the  thermoviscoplastic  material  considered  in 
[  2  ]  .  The  material  model  utilized  here  is  applicable  for  use  with  a  wide 
variety  of  crystalline  metals  at  temperatures  above  half  their  melting  point 

In  order  to  make  the  above  equations  more  tractable,  we  consider  the  case 
of  a  homogeneous  isotropic  thermoviscoplastic  uniaxial  rod.  We  postulate  that 

°ij  -  on  -  a  (1.7) 

*ij  ”  €n  “  * 

A jkl  “  ^llll  *  E 

*ij  “  alij  “  al 

“ij  *  «11  "  « 

*Tij  -  -  /  -  oiT  -  TR ) 

We  shall  also  assume  for  simplicity  that  E ,  a,  k  ,  and  p  are  constant. 

Substituting  the  relations  (1.7)  into  (1.2)  we  have  the  simplified  equation 

£(€  —  <*!  +  or T R  )«!  +  E  oP’TT  —  E  aT  t  —  pcv  T  +  k  A T  -  0 .  ( 1.8) 

The  stress  -  strain  relation  reduces  to 

o  -  £U-  aj  -otiT  -TR)l  (1.9) 

We  also  assume  infinitesimal  deformations 

€  -  ux  . 


We  therefore  have  the  uniaxial  model 
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ax  -  0  (1.10) 

E(.€  —  +  aT  g  )ciq  +  E  oP'TT  —  E  aT  e  —  pcy  T  +  k  bT  *  0  (1.1 1) 

a,  -  fiio(j,T,d  (1.12) 

a  -  £[€  —  «!—  o&T  —  r*)]  (1.13) 

where  the  infinitesimal  strain  is  defined  by 


In  this  simplified  one-dimensional  model,  all  variables  (except  for  T )  are 
independent  of  r  and  9. 

Given  initial  conditions  for  at  ,  T  ,  and  u,  as  well  as  boundary  conditions 
for  T  and  u ,  this  is  a  we  11- posed  mathematical  problem.  It  has  been  shown  in 
[  2  ]  to  be  thermodynamically  consistent  as  well. 

We  make  the  following  observations  at  this  time.  First,  (1.10)  implies 
that  o  is  a  function  of  time  alone.  Secondly,  in  the  special  case  where  the  f , 
do  not  explicitly  depend  on  F  or  €  but  only  on  (i  ,  a)  ,  then 
ocj  -  oi j{t  ;  o(t )).  If  a  were  known,  then  (1.8)  becomes  a  single,  non-linear 
parabolic  differential  equation  for  the  temperature  T  in  terms  of  a.  This 
analysis  has  been  carried  out  in  l  3  ] .  Following  this  approach,  we  will  obtain 
the  asymptotic  behavior  of  T  directly  from  the  system  (1.10)-(1.13)  . 

Section  2.  Radiative  and  convective  boundary  conditions. 

Using  (1.9)  we  can  simplify  (1.1 1)  to 

oax  —  adT  —  pcvT  +  k  AT  -  0  (21) 

In  order  to  include  the  thermal  radiation  boundary  conditions  along  the  lateral 
surface  of  the  rod,  we  integrate  (21)  over  the  cross-sectional  area  of  the  rod, 
which  we  shall  assume  to  be  constant.  Letting  x  denote  the  axial  coordinate, 
we  have 

A  oax  —  AabT  —  A  pcvT  +  kATxx  - 
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where  qn  is  the  normal  component  of  the  thermal  flux  along  the  lateral  sur¬ 
face  of  the  rod,  C  is  the  circumference  of  the  rod,  and  A  is  the  cross-sectional 
area  of  the  rod.  7  is  now  interpreted  as  a  radial  average  of  the  temperature 
over  the  cross-section  of  the  rod. 

It  is  clear  the  term  Cqn  acts  as  a  “sink"  for  thermal  energy.  In  deep  space, 
we  model  the  thermal  flux  due  to  radiation  by  the  Stefan  -  Boltzmann  radia¬ 
tion  condition 

qn  -  -  kse{T*-T$)  -  kseT*-Q  (2.2) 

on 

where  TD  -  0°A'  is  the  deep  space  ambient  temperature,  and 
ks  -  5.775X10-11  is  Boltzmann’s  constant  e  is  an  order  one  constant,  called 
the  emissivity,  which  measures  the  effect  o'"  the  surface  coating  on  black  body 
radiation,  and  Q  measures  the  solar,  earth,  and  deep  space  thermal  flux 
incident  on  the  rod. 

We  linearize  (2.2)  about  7  -  T x  to  obtain  the  convective  boundary  con¬ 
dition 

qn  -  kseTg4  +  0(7  -  7* )  -  Q  (2.3) 

with 

&~~4kseT/(  . 

By  convexity, 

kseT*  >  ks  eTf  +  4 k,  eT^T  -  7 R  )  (2.4) 

if  7  >  7*.  If  we  therefore  choose  3  -  4kseT £  in  (23),  less  heat  is  radiated 
away,  and  we  obtain  an  upper  bound  for  the  temperature  increase.  On  the 
other  hand,  if  we  choose  0  large  enough  so  that 


£ 


.V 


kseT £  +  0(7  —7*)  >  kseT* 


(2.5) 
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over  the  range  of  temperatures  realized,  then  we  obtain  a  lower  bound.  The 
sharpest  lower  bound  of  this  form  is  obtained  when 

ks  eTj  +  {3LT  -  Tr  )  -  ks  eTA 

at  T  “  T  These  observations  follow  directly  from  the  maximum  principle. 
Section  3.  Further  reduction  of  the  equations. 

From  the  boundary  conditions  (2.3),  introduced  in  the  previous  section. 


we  have 


-pcvT  +  kTxx  —  abT  +  oa^  -  +  ~^qn 

A 


+  £  kseTj?  +  0(r  —Tr)  —  Q 


(  For  a  solid  rod,  R  -  radius  of  the  rod. )  We  write  this  as 


pcvT  -  kTxx  +  abT  +  fg.iT  ~  TR  )  -  (3.2) 

A 

o«i  -  lkseTR4  +  W 

Introducing  a  new  dependent  variable  Oix  4  )  =  exp(— —  TR)  , 

pcv 


we  have 


where 


*  -  -L*  +  »  e  -  Fb) 

9*  pcy  ax2  pcv/? 


Fb)  -  expti^  ^1^12 +  ^Q_ 
pcv  pcv  R  pcv 


R  pcy 


■eT R  -  a°-TR 


The  modified  temperature  difference  9  has  been  introduced  previously  in  [  3  ] 
.  Since 


we  have  that  0  *'■*  T  —  Tp.  If  we  integrate  both  sides  of  (3.3)  with  respect  to 
t,  and  divide  by  t,  we  obtain 


t  “KOGc  ,t  )  —  0(x  ,0))  -  — —  t  -1[  f  0{x  ,r)d  r]x 

(x=v  o 

t  t 

+  [9(x  ,r)rfr  -  t~1[F(r)dT 

OCv  R 


As  r— *oo  ,  with  0  bounded,  this  reduces  to 


where 


R  pCy 


<  e  >  -  <&> 


<  F  > 


<  0  >  (x  )  -  lim  t~l  [9(x  ,r)d  r 

l  — »CO  Q 


<  F  >  -  llmr-/«p(«2<d)  **±1  +  22_  (3.6) 

t  — *oo  q  pCy  pCy  A  pCy 

^^5  4  CKO  m  » 

—  - eT  p  — - /  o  cfr 

PCy  pCy 

-  ltor-/exp(^)  g(±ltr)  -  »!_*•/  +  _S_Lr 

f  —*00  ^  pCy  pcv  pcv  R  P^v 

are  the  asymptotic  mean  values  of  0  and  F  respectively. 

Note  that  as  R— >0,  <  F  >  becomes  negative,  and  therefore  T  <Tp. 
Conversely,  as  R  — >00,  <  F  >  becomes  positive,  and  therefore  T  >  Tp.  Intui¬ 
tively,  thin  rods  radiate  away  energy  through  the  lateral  surface  area  more 
quickly  than  thicker  rods. 
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In  the  materials  under  consideration,  oa1  >  0  .  Since  the  integrands  of 
(3.5)  ,  (3.6)  are  bounded,  the  above  limits  exist-  For  simplicity,  consider  the 
following  boundary  conditions  for  the  temperature  T ,  TiOj )  -  TR  and 
Tx  (l,f )  -  0.  The  second  boundary  condition  results  from  symmetry,  for  a 
bar  of  length  two  units.  (More  general  boundary  conditions  will  be  considered 
in  the  next  section.)  This  implies  that  0  satisfies  the  boundary  conditions 

<  0  >  (0)  -  a  (3.7) 

and 

<  0  >  '(1)  =  0  (3.8) 


The  solution  of  (33),  subject  to  the  boundary  conditions  (3.7),(3«8)  is 


<  0  >  (x ) 


pcvR  cosh\(l— x) 

-VT<F>  1  ~  coshX 


(3.9) 


where 


X 


r  2/3  m/2 
Rk 


Note  that  if  we  ignore  the  spatial  variation  of  the  temperature,  that  is 
assuming  that  Txx  -  0,  we  have 

.  ,  .  pcv R 

<e>(x)  -  -^-<F  > 

which  is  a  good  approximation  to  (3.9)  away  from  the  ends  of  the  rod,  if  k  is 
large.  This  accounts  for  the  close  approximation  observed  between  the  spa¬ 
tially  homogeneous  solution  and  the  spatially  varying  solution  with  boundary 
layer  in  [  4  ] .  We  remark  that  (3.9)  yields  x  ~X~1  as  a  measure  of  the  asymp¬ 
totic  thickness  of  the  boundary  layer  for  the  boundary  conditions  considered 
above. 

Since  p,Cy  ,  1 3  ,  R  ,  and  k  are  all  known  ,  the  spatial  dependence  of 
<  0  >  (x)  is  determined.  Only  the  magnitude  is  a  function  of  the  loading 
history,  through  <  F  >  .  The  fine  structure  of  the  material  properties  has 


been  averaged  out,  and  the  asymptotic  form  of  the  history  dependence  takes 
the  form  of  a  single  scalar  quantity  <  F  >  .  Many  internal  state  variable 
models  for  the  stress  history  can  therefore  be  considered.  The  distinctions 
between  the  models  are  lost,  or  more  precisely  averaged  out,  in  this  analysis. 
The  crucial  term  to  estimate  therefore  is  <  F  >  .  If  ait )  is  periodic,  and  if 
<*!  is  a  periodic,  then  (3.6)  reduces  to 


<  F  > 


r 

'l  fexpt 


oeo  ^  _  2 

pCy  pCy  R  pCy 


.kseTBA\dr  (3.10) 


where  P  is  the  period  of  o  and  ij. 

A  case  of  particular  interest  is  the  Bodner  model  [  5  ]  for  metals  ,  in 
which 


or,  -  *  ZJo^exp  -  As^-I— |2n 

1  73  jof  2n  1  a  1 


(3.11) 


At  saturation,  the  drag  stress  a2  -  a2  is  constant,  and  <  F  >  becomes  a 
functional  only  of  a. 


<  F  >  -  P-'fex p 

J0  pcv 


(3.12) 


^ £>oKr)|exp  -  A^-|-T^r|2n  “  — M?*4  +  ,  d  r 

/3  pcv  1  1  2 n  'TtrT  R  pcy  *  Rpcy 


Since 


<0>  -  lim  t  ~l  fOix  ,r)d  t  «  limf  1  fexpt-^^-Xr  —TR)d  t 

<-  -  0  0  PCV 


we  can  recover  a  mean  temperature  T  by 


T  -  <  9  >  (jc  )  limr  1  f expaa^T^d  t  +Tt 

1-00  J0  pCy 


(3.13) 


Fixing  the  geometric  parameters  {  C  ,A  },  and  fixing  the  material  param¬ 
eters  {  a,p,cv  ,D0,n  ,&,k  },  one  can  compute  the  dependence  of  <  0  > 
and  T  on  o(t)  and  Q  .  With  a  periodic,  (3.13)  simplifies  as  before  to  an 
integral  over  the  period.  A  comparison  of  analytical  results  with  experimental 
and  numerical  data  is  contained  in  section  7  for  a  particular  case  of  periodic  o 
,  as  illustrated  in  fig.  2. 

Because  (3.3)  is  linear  in  0,  we  can  obtain  explicit  representations  for  0  in 
terms  of  series.  Information  on  the  growth  of  boundary  layers,  and  an  explicit 
compliance  relationship  between  strain  and  stress  will  also  follow  from  an 
analysis  of  (3.3).  First  we  derive  a  series  representation  for  0. 

Section  4.  Series  expansion  of  the  modified  temperature  difference. 

We  consider  the  solution  of  problem  (3.3) 

+  J»  t-Fbi 

df  pCy  hx2  PCyR 


with  initial  conditions 


0(x  ,0)  -  0 

and  slightly  more  general  boundary  conditions  (  with  a  parameter  /3  ) 

-  J50  -  0 


at  x-0  and 


at  x-1.  We  write  the  solution  of  this  problem  in  the  form 


0U,t)-  £0n(r)*nU). 

n  =1 


Substituting  this  into  (33),  we  have 


I  QnXt)*n{x)-  -F(f)( 4.1) 

n  =1  PCV  P^V  * 


with  boundary  conditions 


xe^o^ui-o 

n  =1 


£0n(/)^n'(l)-O 


£  e„  (f )  k  t/fn  '(o)  -  (o)  I  -  a 

n  =1  1  1 


We  solve  (4.1M4.4)  by  assuming  that 

fl„  Xi )  +  (< )  +  -  J0  e„  (< )  -  dn  F  (r ) 

PCy  pC y  K 


\ K  "(x )  +  ca*\l}n  (x )  -  0 


Therefore, 


£  en  •(/  0c )  -  JLs.  (r  ••(* )  +  -I®  en  (<  c* ) 

n  =1  PCV  PCV/C 


This  will  satisfy  (4.1)  if  and  only  if 

ZdnMfn  (x)-l  (4.7) 

n  =1 

on  the  interval  0  <  x  <1.  Equations  (4.3)  ,  (4.4)  ,  and  (4.6)  determine  co„ 
uniquely.  One  can  easily  show 

*l>„(x )  -  cos(l— x  ^ 
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Asymptotically , 


r <*n  -  cotton 


ton  *n  it  , 


We  note  that  ton  is  nonzero,  so  that  \Utx )  -  1  is  not  an  eigenfunction. 
Therefore,  (4.7)  determines  the  coefficients  dn  by 

£rfncos(l -*)<«>„  -  1  (4.9) 

n  =1 

The  eigenfunctions  xl/„(.x )  «  cos(l—  x  ]o„  area  complete  orthogonal  basis,  con¬ 
sequently 


4sincon 

2ton  +  sin2o>„ 


(4.10) 


Asymptotically  , 


From  (4.2) 


2sino>„  23 
<Ai  kn2ir2 


G„(0)-0. 


(4.11) 


The  solution  of  (3.3)  can  be  written  down  concisely  as 


0(x ,  t )  -  £  dn  j*exp(— y2{t  —s  ))F (s  )ds  Los(l— x  )o>n  (4.12) 


n  =1  I  0 


where 


mu 


and  dn  satisfies  (4.1  OX 

Section  5.  Boundary  Layer  Model. 

In  order  to  examine  the  behavior  of  the  boundary  layer,  we  expand  the 

function  ,  t )  in  terms  of  the  small  parameter  y2  -  — —  .  We  define  the 

pcv 

stretched  coordinate  £  -  y~lx .  This  implies  that  dx  -  y-I9f,  and  that  there¬ 
fore 

e,  e-/^)  (5.1) 

pcv  R 

e(y£.0)-0 

*y'1e£(0,r)- 3e(0,r  )-0 
In  the  usual  way,  we  write 

«{■<)-  £yM{,») 

n  =0 

For  n-0  this  reduces  to 

%  1  ~  %  u  +  — •Fit)  (5.2) 

pCy  K 

e0(  • .  0)  -  0 


with  boundary  conditions 


90,{(oo,  t )  -  0 


e0.*(o,*)-o 


which  implies  that  0O  -  0O(? ),  and  that 


0o-eo(O-  /F(r)exp(— i*L(f 
o  Pcv  R 


-  r))d  r 


Equation  (5.3)  describes  the  evolution  of  the  modified  temperature  difference  0 
under  insulated  boundary  conditions  at  the  ends  of  the  rod.  For  n>l  ,  we 
obtain  a  coupled  set  of  linear  differential  equations 


Qn.t  ~  ®n.tt  "  0 


0„(.  ,0)«0 


*)  -  0 

*0B,*(aO- J59n_i(ftf)-0 


This  can  be  solved  using  Laplace  transforms  ,  which  yields 


*<«•'>-  “lv?rexp("lr)‘e"-l(a') 


Note  that  if  Q0(t )  >  0  then  each  term  alternates  in  sign. 


The  first  order  approximation  is  therefore 


0(x  ,t )  0o(f  )  +  y0i(jc  ,t ) 


_  t  2 

0O(*  )  ~  ~  f~r — exP(  “  — — )0o(f  ~  f)d 

k  %  yJVT  Ay  T 


The  integrand  is  negligible  except  when  is  of  order  one.  This  leads  to  a 

y  t 


boundary  layer  which  is  described  for  small  time  by 


~y  \/T  . 


We  also  obtain  the  following  result 


U  )-e0(t )  -  y  If  ^i-GoC  t  —r)d 


which  to  leading  order  computes  the  temperature  as  a  function  of  time  of  the 
left  endpoint  of  the  bar.  The  first  order  correction  is  in  fact  monotonic  in  x  if 
90{t )  >  0  .  Equation  (5.7)  shows  that  the  first  order  term  overcorrects  and 
that  the  next  term  (of  opposite  sign)  attempts  to  compensate  for  this. 

This  singular  perturbation  analysis  predicts  the  behavior  of  the  rod  for 
small  time,  in  contrast  to  the  analysis  of  section  3  which  predicts  the  long¬ 
term  average  behavior. 

Section  6.  Integral  representation  of  9{x ,  t )  by  means  of  a  Green’s 
function.  Derivation  of  the  Compliance  Relation. 

By  taking  the  Laplace  transform  of  (3.4)  with  the  boundary  conditions 
introduced  in  the  previous  section,  we  obtain  the  following  expression  for  the 
Laplace  transform  of  9 


(Kx ,  s  )  -  (s  +  —?.--)  1F  (s  )G  (x  ,  s  ) . 
pcy  R 


0OU ,  s  )G  (x  ,  s  ) 


where  90  was  defined  in  (5.3)  and 


G  (x ,  s )  «  1  — 


)3cosh 


j(  1  — X  )\/o> 


k  vcosinhvca  —  /3coshvc o 


,  to -is  +  J^X^L)  (6.2) 


pCy  R 


Inverting  this,  we  obtain 


@(x  JL )  -  G  (x  4  )  *  90(t ). 


Using  (1.9)  and  the  definition  of  9  we  have 


€  -  £  lo  -I-  <*i(f  ;  o(t ))  +  a  exp(—  )9{x  jt  X 

pcy 


V  ,* 
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In  the  simplest  case,  oq  is  obtained  by  integrating 

a,  -  /(t;o(f)X 

and  9  depends  on  a  through  F .  Equation  (63)  yields  the  compliance  relation¬ 
ship  between  the  strain  €  and  the  stress  o.  Further  details  may  be  found  in  [ 
3]. 

Section  7.  Numerical  Results 

In  order  to  verify  the  accuracy  of  the  asymptotic  result  (3.9)  ,  we  com¬ 
puted  upper  and  lower  bounds  for  the  asymptotic  average  modified  tempera¬ 
ture  difference  for  the  parameters  given  in  [  6  ]  with  the  results  shown  in  fig. 
1.  The  value  for  the  emissivity  is  0.85  and  -  346.8  MPa  corresponding  to 
a  1%  maximum  deformation  at  5  hz. 

The  cyclically  saturated  stress-strain  curve  used  to  compute  <  F  >  for 
the  case  -  3363  MPa  is  given  in  fig.  2.  The  values  for  omjLX  ,  a2  ,  and  the 
remaining  parameters  needed  to  compute  <  F  >  and  <  9  >,  were  taken 
from  numerical  data  supplied  from  reference  [  6  ] . 

The  parameters  in  [  6  ]  describe  a  hollow  cylindrical  bar,  of  uniform 
cross-  sectional  area  equal  to  1  in.2  ,  orbiting  at  an  altitude  of  4080  km., 
painted  with  a  high  emissivity  coating  (IITRE-S13GLO)  white  paint,  with  full 
exposure  to  the  sun.  Under  zero  applied  stress,  equilibrium  of  thermal  flux 
occurs  at  72  deg.  Fahrenheit  (295  deg.  Kelvin)  for  which  the  thermodynamic 
parameters  of  the  metal  6061-T6  aluminum  were  experimentally  determined  [ 
7). 

We  also  can  compute  the  behavior  of  the  thermal  boundary  layers  near 
the  ends  of  the  rod.  In  equation  (3.3)  the  diffusion  constant 

y  -  (_^_)1/2  -  23x1 0-4  .  This  quantity  typically  is  of  the  order  of  10-3  to 
pcv 

10~3  in  dimensionless  units.  This  leads  to  very  sharp  boundary  layers,  and 
exponential  decay  to  the  asymptotic  temperature  state. 

We  also  mention  that  other  constitutive  laws  than  (3.11)  have  been 
derived  which  describe  the  behavior  of  metal-matrix  materials.  Since  the 
parameters  are  experimentally  determined  by  a  parameter  fitting  procedure. 
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essentially  the  same  results  ■will  be  obtained  for  the  asymptotic  temperature 
rise. 

The  interested  reader  is  referred  to  [  8  ]  for  a  discussion  of  several  of  these 
models. 

Section  8.  Wave  Dispersion  in  the  Isothermal  Case. 

Consider  the  equations  of  motion  for  an  uniaxial  rod  of  infinite  length 
with  constant  cross-section,  temperature,  and  density  : 


where 


PUtt 


3a 


3x 

a  -  E(e  —  aj) 
3al  r(  \ 


€  - 


3u 

dx 


(8.1) 

(&2) 

(a3) 


is  the  infinitesimal  strain,  and  aj  is  the  internal  state  variable  describing  ine¬ 
lastic  strain.  These  equations  are  limiting  cases  of  equations  considered  in  sec¬ 
tions  1  -  6  ,  with  inertial  effects  added. 

We  can  combine  equations  (8.1)  and  (8.2)  into 

putt  -  ±[E  (|i  -  a,)]  -Eu^-E  ou  (8.4) 

Writing  (83)  and  (8.4)  as  a  system  of  3  first  order  partial  differential  equa¬ 
tions,  we  let  ux  -  Uj,  u2  -  and  u3  -  orj.  Consequently 

pu2't  -  Euhx  -Euyx  -  £(u,  -  u3)_  x  ( 8.5 ) 

u3  t  -  f{E{ux  —  u3))  (8.6) 

with  the  compatibility  condition 

“u  ““2,j  (8.7) 


La 


ii 


v 


N  , 


..Xk  .-J 


Since  the  quantity  Uj  —  u3  figures  predominately  in  equations  (8.5)  and 

(8.6),  we  define  new  independent  variables  :  w  j  -  u  t 

and  we  have 

-  Up  m>2  -  U2  ,  W3  -  u3 

pw  lt  -  £w  u 

(8.8) 

-  f  (E^i) 

(8.9) 

-  w2  :c 

Eliminating  w  3il  in  (8.10),  by  (8.9),  we  obtain 

pw2fl  -  £wu 

™3j  “  f(Ewx) 

-  w2jc  -  -  f(Ewx) 

which  we  write  in  canonical  form 

(8.10) 

vu  “  w2a  -  “  f(Ew ,) 

(8.11) 

h 

1 

'o  |  th 

Vi 

1 

o 

(8.12) 

«  +  f  {Em-  ,) 

(8.13) 

The  characteristics  of  (8.1 1M8.13)  are  given  by 


dt 

r 

where  A.  =  0,  ±(_),/2  .  We  therefore  have  3  constant  characteristic  wave 
P 

speeds.  As  we  shall  see,  because  of  the  nonlinearity  induced  by  ax  ,  waves 
propagate  with  a  continuum  of  wave  speeds. 

We  seek  traveling  wave  solutions  of  (8.1 1M8.13)  ,  that  is  we  seek  solu¬ 
tions  of  the  form  v,  -  M-t{x  —  ct )  -  w,  (£)  where  c  -  constant  is  the  wave 
speed.  Substituting  this  in  the  system  of  partial  differential  equations,  we 
arrive  at  a  first  order  system  of  nonlinear  ordinary  differential  equations. 


T 


-w3-  +/(£wj)  (8.16) 

where  dot  signifies  differentiation  with  respect  to  £  »  x  —  ct .  Equations 
(8.14M8.16)  decouple  very  nicely  by  using  (8.15)  to  eliminate  w2  from  (8.14). 
We  therefore  obtain 

—  cvj  +  -  — /(Ew,) 


or 


(c2  —  — )w  i  =  c/  (Ew  j) 
P 


which  reduces  to 


(c2- -£XEwj)»  Ec/(Ew!)  (8.17) 

P 

The  quantity  Ew  j  is  of  particular  interest,  since 

Ew  j  -  Etu  !  —  u3)  -  E(ux  —  aq)  =  a 
Consequently,  (8.17)  can  be  written  as 

^  °  c£  „  /  (a)  (8.18) 

c2-^ 

P 

2  E 

If  c  =  —  then  the  wave  is  traveling  at  sonic  speed  and  b  may  be  unbounded. 
We  consider  only  the  subsonic  case  in  this  paper. 

It  is  interesting  to  note  that  the  equation  describing  the  propagation  of  the 
stress  wave  a(x  —ct )  decouples  from  the  rest  of  the  equations.  Equation  (8.18) 
may  be  solved,  for  a  given  f ,  for  a  continuum  of  wave  speeds  c.  We  inter¬ 
pret  the  quantity 


as  a  nonlinear  amplitude  .  In  the  Bodner  model,  for  example, 

(8-19) 

where  Z)0,  n ,  and  a2  are  material  parameters. 

Note  that  the  f,  given  by  (8.19),  has  the  property  that 
oat  =  a  f  (o)  >  0.  This  implies  that  any  traveling  wave  solution  of  (8.18) 
o(x  —ct )  with  compact  support  must  be  discontinuous.  This  follows  by  noting 
from  (8.18)  that 

y2J*o2-o^L- — SEl—o  f  (a)  ~  A{c)o  f  (a)  (8.20) 

dS  c2_£. 

P 

For  a  given  c ,  A  (c  )  is  either  positive  or  negative  but  has  fixed  sign.  Conse¬ 
quently  the  amplitude  of  a  is  monotonically  increasing  or  decreasing,  which 
implies  that  a  continuous  stress  field  a  may  vanish  at  most  once.  Given  a 
smooth  initial  stress  field  with  compact  support,  a  single  wave  speed  is  not  pos¬ 
sible,  and  the  disturbance  propagates  in  a  continuum  of  admissible  wave  speeds 
with  nonlinear  interactions. 

The  inequality 

aafi  -  of  (a)  >  0 

which  is  responsible  for  the  dispersive  effects,  is  also  responsible  for  the 
conversion  of  strain  energy  into  heat  in  thermoviscoplastic  materials  (cf.  3.9  , 
3.10)  .  This  dispersion  also  tends  to  drive  an  initial  impulse  towards  a  quasi¬ 
static  state  where 


pu[t  -0,-0 


justifying  the  assumption  of  spatial  homogeneity  in  one  dimension. 


Note  that  we  can  recover  the  initial  displacement  field  by 

u0'(x )  - +  otjCj:  X)).  (8.21) 

For  t  >  0  ,  we  can  recover  u  (x  ,t )  by  integrating 

a  . 

“x  -  -j-  +  <*i  • 

Note  also  that  the  wave  speed  c  -  0  is  prohibited,  for  if  c  -  0  then  a  -  a(x ) 

which  implies  that  -  0  which  can  only  happen  if  /  -  0  identically  in 

at 

(8.31  Therefore  for  smooth  traveling  wave  solutions  of  (8.1)  -  (8.3) 

0<|c|  <(*)">  (8.22) 

P 

Section  9.  Conclusions. 

In  terms  of  the  modified  temperature  difference  9  ,  we  have  derived 
upper  and  lower  bounds  on  the  asymptotic  average  temperature  difference. 
With  (2.4)  we  obtain  upper  bounds,  and  with  (2.5)  we  obtain  lower  bounds 
via  the  maximum  principle.  The  lower  bounds  are  sharper,  and  with  the 
parameters  described  in  section  7  ,  a  frequency  of  5  hz  (corresponding  to  omax 
of  346.8  MPa)  an  upper  bound  of  843°  K  and  a  lower  bound  of  61.9°  K  on 
the  temperature  change  are  obtained.  The  upper  and  lower  bounds  obtained  in 
this  paper  are  in  close  agreement  with  the  numerical  results  obtained  in  refs.  [ 
4  ]  and  [  5  1  Temperature  increases  on  this  order  can  lead  to  material  failure 
of  AL  6061-T6,  suggesting  that  material  inelasticity  be  used  with  caution 
when  operating  near  yield  strength  for  metal  matrix  materials. 

Explicit  series  representations  for  0(x  jt )  in  terms  of  the  stress,  through 
Fit ),  are  derived.  The  compliance  relation  involves  convolutions  as  expected  , 
but  is  not  a  pure  convolution  due  to  the  presence  of  the  internal  state  variable 
describing  inelastic  strain. 

The  asymptotic  average  modified  temperature  difference  ,  given  by  (3.9) 
has  a  certain  canonical  form.  The  amplitude  depends  on  geometric  and 


material  parameters  (constants)  in  a  simple  way,  and  depends  on  the  particular 
model  used  only  through  the  scalar  term  <  F  >  .  The  longitudinal  varia¬ 
tion  is  nearly  constant,  except  near  the  ends  of  the  rod  where  boundary  layers 
form.  These  boundary  layers  are  described  in  Section  5. 

Boundary  layers  are  expected  to  form  in  the  radial  direction,  but  these 
will  serve  primarily  to  reduce  the  temperature  on  the  outside,  leading  to  less 
of  a  heat  loss  than  otherwise  predicted.  Therefore  the  temperature  differences 
in  the  interior  will  be  larger  than  the  average  temperature  difference  in  the 
cross  section.  Boundary  conditions  at  the  ends  of  the  rod  do  not  substantially 
affect  the  thermomechanical  response  outside  the  boundary  layer  for  the  struc¬ 
tural  configurations  considered  in  this  paper. 

The  uniaxial  (one-dimensional)  case  we  have  considered  in  this  paper  is 
much  simpler  than  the  full  three-dimensional  problem  (1.1)  -  (1.6).  This 
simplification  allows  us  to  obtain  much  stronger  results  however.  We  incor¬ 
porate  time-dependent  effects,  full  thermal  -  mechanical  coupling,  and  non¬ 
linear  black  body  radiation  conditions  to  obtain  rigorous  bounds  on  the  asymp¬ 
totic  mean  temperature  rise.  In  those  cases  where  longitudinal  effects  dom¬ 
inate  axial  effects  and  the  resultant  quantities  are  approximately  uniaxial,  the 
analysis  presented  in  this  paper  is  expected  to  be  qualitatively  correct. 

Finally,  the  dispersive  effects  due  to  material  inelasticity  are  discussed  in 
Section  8.  It  is  found  that  dissipation  and  dispersion  effects  generally  inhibit 
the  formation  of  traveling  waves,  and  tend  to  drive  the  impulse  into  a  quasi¬ 
static  state,  thereby  justifying  the  neglect  of  inertial  effects  in  sections  1-7. 
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Abstract 

A  finite  element  model  Is  outlined  for  an 
integrated  thermoelastlo  analysis  of  large 
composite  space  structures.  The  model  allows  tor 
temperature  gradients  within  structural  member 
cross-sections  and  for  bending  of  the  members 
themselves.  Nonlinear  effects  such  as  radiation 
boundary  conditions  and  temperature  dependent 
material  properties  are  also  Included.  Once  the 
model  Is  outlined,  a  preliminary  Investigation 
into  the  importance  of  thermally  Induced  forces 
and  moments  Is  carried  out.  The  problem  chosen 
Is  that  of  a  long  cantilevered  lattice  beam  In  a 
geosychronous  orbit.  For  the  structure  and 
loading  chosen,  no  significant  dynamic  responses, 
such  as  vibration,  occurred.  In  addition, 
thermally  Induced  axial  forces  were  the 
predominant  type  of  loading.  For  this  problem, 
thermally  Induced  moments  could  be  neglected. 
The  magnitude  of  axial  stresses  generated  by  the 
transition  from  shadow  to  sunlight  Is  on  the 
order  of  30J  of  yield  stress. 


Nomenclature 

A  cross-sectional  area  of  a  structural  member 

B  boundary  of  cross-seotlonal  area 

cp  specific  heat 

E  Young's  Modulus 

G  modulus  of  rigidity 

I  mass  Inertia 

o 

Iy,Is  moments  of  Inertia  about  y  and  z  axes 

J  polar  moment  of  inertia 

k  thermal  conductivity 

L  length  of  a  structural  member 

T 

M  thermally  Induced  moment 

n  direction  normals 

T 

P  thermally  Induced  axial  force 

q  normal  Incident  flux 

T  temperature 

T0  temperature  In  an  unstrained  state 

Tp  radiation  reference  temperature 

t  time 

u,v,w  displacements  In  coordinate  directions 
x.y.z  coordinate  directions 

*Graduate  Resear oh  Assistant,  Aerospace 
Engineering,  Student  Member  AIAA 
••Associate  Professor,  Aerospace  Engineering, 
Member  AIAA 

•••Professor  and  Head,  Aerospace  Engineering, 
Associate  Fellow,  AIAA 


coefficient  of  thermal  expansion 

e  emlsslvlty 

3  rotation  about  x  axis 

p  density 

a  Boltzman's  constant 

Introduction 

In  the  next  few  decades  large  space 
structures  will  be  placed  Into  earth  orbit. 
Because  these  vehicles  will  not  face  the  launch 
environment  but  are  to  be  either  deployed  or 
constructed  from  materials  carried  into  orbit, 
design  will  be  based  on  criteria  previously 
considered  secondary.  These  criteria  combined 
with  new  performance  requirements  will  result  in 
large  lattice-type  structures  which  utilize  high- 
strength,  low-density  materials  such  as  graphite 
flber/polymerlc  matrix  composites. 

Under  anticipated  thermal  loading 
conditions,  structural  members  made  from  advanced 
fiber /matrix  materials  respond  quite  differently 
compared  to  those  made  from  the  more  traditional 
metallic  materials.  Past  research  indicates  that 
a  graphite/epoxy  structural  member  modeled  as  a 
slender  thln-walled  tube  experiences  a 
significant  temperature  gradient  around  Its 
perimeter.1  In  addition,  such  a  member  will  have 
a  negligible  temperature  gradient  along  its 
length.1  Both  these  responses  are  due  to  the  low 
thermal  conductivity  found  In  fiber/matrix 
materials. 

A  large  temperature  gradient  through  a 
structural  member's  cross-section  produces 
bending.  This  Is  Important  for  two  reasons. 
First,  bending  reduces  the  maximum  allowable  load 
a  member  can  sustain.  Second,  bending  may  lead 
to  fatigue  which  is  important  In  predicting  the 
long  term  behavior  of  the  material.’ 

In  this  paper,  the  development  of  a  finite 
element  model  capable  of  performing  an  Integrated 
thermoelaatic  analysis  Is  outlined.  This  model 
Is  specifically  designed  for  lattice-type 
structures  made  from  low  conductivity  materials, 
wherein  the  temperature  distribution  within  a 
structural  member  varies  through  Its  cross- 
section  but  not  along  its  length.  In  addition, 
preliminary  studies  are  carried  out  to  determine 
the  significance  of  thermally  Induced  bending  and 
extension.  For  this  purpose  a  typical  space 
structure  is  developed  and  Its  thermal/structural 
responses  examined  for  two  different  load  cases. 
Further  details  are  contained  In  reference  4. 

Background 


A  general  thermoelastlo  analysis  or  a  large 
space  structure  is  complicated  by  several 
factors.  First,  there  is  the  coupling  between  the 
temperature  field  and  the  displacement  field. 
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boundary  conditions  and  temperature  dependent 
material  properties.  Third,  the  thermal  loading 
is  constantly  changing  due  to  varying  earth- 
structure-sun  orientation.  Finally,  geometrical 
factors  such  as  shadowing  and  lnterlement 
radiation  can  exist. 

In  the  model  developed  herein,  several 
simplifying  assumptions  are  made.  To  begin  with, 
the  temperatures  are  assumed  to  be  Independent  of 
''^formation,  that  1s,  uncoupled.  This  allows  the 
temperature  field  to  be  solved  for  first,  then 
used  as  Input  to  the  structural  analysis. 
Second,  all  structures  are  assumed  to  be  of  open 
lattice-work  construction  with  long  thin  members 
sparsely  located.  This  allows  for  the  omission 
of  the  lnterlement  radiation  and  shadowing.  The 
structural  members  are  modeled  as  thln-walled 
cylinders  of  constant  cross-section  and  material 
properties.  Furthermore,  because  temperatures 
within  these  members  vary  through  their  cross- 
section  and  not  along  their  length,  each  member 
may  be  treated  as  an  Isolated  body,  absorbing 
thermal  radiation  and  in  turn  emitting  Its  own 
radiation.  See  Fig.  1. 
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dependent  material  properties,  namely  k  ,  k  , 
and  c  .  Solving  the  nonlinear  heat  t Pans fir 
problem  in  a  single  cross-section  per  member  is 
prohibitively  expensive.  However,  it  Is  expected 
that  the  repetitive  nature  of  the  proposed 
lattice  structures  will  allow  a  greatly  reduced 
thermal  analysis  wherein  the  results  from  a  few 
selected  cross-sections  are  used  throughout  the 
structure. 


STRUCTURAL 

MEMBER 


Figure  1 .  A  Structural  Element  in  Space 
Environment 

The  thermal  loading  Is  also  simplified  by 
requiring  all  structures  to  be  In  a  high  earth 
orbit.  At  high  altitudes  earth  emitted  and  earth 
reflected  radiation  are  negligible.  Furthermore, 
If  structures  are  assumed  to  have  a  space-fixed 
orientation,  the  solar  flux  Is  both  constant  In 
magnitude  and  direction  during  the  sunlit  portion 
of  an  orbit. 

Model  Development 

The  finite  element  model  oonslsts  of  five 
sections.  The  algorithm,  shown  schematically  In 
Fig.  2,  Is  as  follows.  On  a  given  time  step,  the 
proper  thermal  loads  are  evaluated.  Then,  finite 
elements  are  used  to  construct  the  temperature 
fields  through  selected  cross-sections.  The  two- 
dimensional  heat  transfer  problem  within  a  cross- 
section  Is  shown  In  Fig.  3  and  governed  by 
equations  (1)  and  (2).  These  equations  neglect 
Internal  heat  generation  but  Include  temperature 


Figure  2.  Algorithm  Schematic 
SOLAR  RADIATION 
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Figure  3.  Heat  Transfer  In  a  Selected  Cross- 
Section 
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Within  a  cross-section  heat  Is  transferred 
by  conduction  around  the  perimeter  and  radiation 
from  the  inner  surfaces.  In  addition,  heat  is 
lost  to  space  through  radiation  from  the  outer 
surfaces.  The  resulting  temperature  fields  are 
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equating  the  flux  and  radiation  boundary  terms  In 
equation  (2)* 


Case  2 

In  this  oase,  a  transition  from  shadow  Into 
sunlight  is  considered.  Here,  the  boom  is 
undeformed  and  at  an  Initial  uniform  temperature 
of  100  degrees  Kelvin.  At  time  t-0,  the  boom 
moves  into  sunlight  and  deforms.  In  a 
geosychronous  orbit,  penumbra  effects  are  Ignored 
and  the  solar  flux  Is  applied  Instantaneously  and 
uniformly  along  the  length  of  the  boom.  The 
solar  flux  Is  considered  to  act  In  the  negative  z 
direction  with  respect  to  the  global  coordinates 
in  Fig.  4.  Its  magnitude  Is  given  as  1.4  kW/m  . 


5,  It  Is  apparent  that  all  twelve  members  can  be 
assigned  to  one  of  tnree  groups  based  on  their 
orientation  with  respect  to  the  direction  of  the 
solar  riux.  This  grouping  Is  shown  in  Table  3. 
It  should  be  noted  that  if  cross-sectional  or 
material  properties  varied  from  member  to  member 
additional  subgroupings  would  be  necessary.  For 
this  simplified  case,  modeling  only  one  cross- 
section  from  each  group  is  needed  for  determining 
the  temperature  distribution  In  a  bay. 
Therefore,  the  thermal  analysis  of  the  entire 
boom  has  been  reduced  to  examining  three  member 
cross-sections. 


Table  3.  Cross-Section  Numbers  for 
First  Bay  Elements 
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Table  2. 

Orbital 

Data 

GEO 

LEO 

Period  (Sec): 

36,400 

5.400 

Altitude  (km): 

35.800 

230 

Umbra  (Sec): 

4,200 

2,200 

Penumbra  (Sec): 

130 

8 

Finite  Element  Model 

Structural  modeling  of  the  boom  Is 
straight-forward.  The  finite  element 
representation  of  the  boom  uses  one  space  frame 
element  per  member.  This  results  In  a  mesh  of  33 
nodes  and  93  elements.  The  mesh  of  the  first  bay 
is  shown  In  Fig.  5. 
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Figure  5.  Finite  Element  Mesh  for  the  First  Bay 


By  comparison,  thermal  modeling  of  the 
structure  appears  quite  complex.  However,  it  can 
be  greatly  simplified.  Because  each  bay  Is 
identical  to  the  next,  and  the  Incident  flux  is 
botn  uniform  along  the  boom's  length  and  constant 
with  time,  only  one  bay  need  be  modeled. 

Thetemperature  distribution  In  the  bay  Is  then 
applicable  throughout  the  boom.  In  addition,  the 
thermal  analysis  of  the  bay  may  be  simplified. 
Upon  examination  of  the  first  bay,  shown  in  Fig. 
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Figure  6.  Finite  Element  Mesh  for  a  Selected 
Cross-Section 


The  finite  element  representation  of  a 
cross-section  Is  shown  In  Fig.  6.  Here  only  half 
the  cross-section  need  be  modeled  if  the  local 
axes  of  the  cross-section  are  aligned  with  the 
solar  flux  vector.  The  mesh  used  consists  of  4d 
constant  thermal  gradient  elements  and  39  nodes. 
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than  converted  Into  moments.  and  axial  forces 
using  the  following  relations  . 
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These  loads  are  Initially  determined  in  the  local 
coordinates  of  their  cross-section.  Later  they 
are  transformed  Into  the  global  coordinates  of 
the  structure  to  be  used  as  Input  for  the 
structural  analysis.  The  struotural  (or 
deformation)  analysis  Is  for  linear  elastic  space 
frame  geometries.  It  is  developed  from 
application  of  standard  finite  element 
formulations  to  the  following  governing  equations 
of  beam  motion7’  . 
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Figure  4.  Ten  Bay  Lattice  Boom 


Table  1 .  Material  Properties  for 

Graphite  AS/Epoxy  Composite 
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Once  the  deformations  are  determined,  time  Is 
Incremented  and  the  process  repeated. 


TOO k  (w/m/K) cp  (J/Kg/K) 


0 

0.0 

0.419 

120 

3.834 

338.0 

(7) 

170 

5.993 

479.0 

220 

8.032 

620.0 

270 

9:T14 

783.0 

330 

10.14 

976.0 

(8) 

400 

11.14 

1080. 

810 

16.98 

1660. 

E  -  4.5  x 

1010  N/m2 

(9) 

G  •  1 .5  x 

1010  M/m2 

p  -  1.633 

x  103  Kg/m3 

a  -  0.91 6 


Problem  Summary 


In  order  to  demonstrate  the  model's  use  and 
as  an  initial  Investigation  Into  the  significance 
of  thermal  loads,  consider  the  following  problem. 
A  long  cantilevered  boom  structure  is  In  a 
geosychronous  orbit  about  the  earth.  During  an 
orbit,  the  boom  moves  from  earth's  shadow  into 
sunlight  and  back  Into  shadow.  While 
transitioning  from  one  thermal  environment  to 
another,  the  boom  undergoes  deformation  and 
possibly  vibration. 

The  structure,  shown  In  Fig.  A,  is  SO  m 
long  and  composed  of  ten  Identical  bays  of 
triangular  cross-section.  Each  bay  Is  5  m  long 
and  2.5  a  on  a  side.  The  structural  members 
themselves  are  thln-walled  cylinders  and,  except 
for  length.  Identical  to  one  another.  The  cross- 
section  diameter  is  0.1  a  and  the  thickness  is 
.004  a.  All  members  are  assumed  to  be  made  from 
Graphite  AS/Epoxy  with  the  properties  given  In 
Table  1.  It  should  be  noted  that  these 
properties  are  for  a  quasi- Iso tropic  lay-up. 
Very  little  Information  is  available  on 
properties  of  complex  lay-ups. 
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Two  load  cases  are  examined.  The  first 
consists  of  the  boom  structure  originating  In 
sunlight  and  moving  Into  shadow.  The  second 
examines  the  same  boom  moving  from  shadow  into 
sunlight. 


For  the  first  case,  the  structure  is 
considered  to  be  In  thermal  equilibrium  and 
stress-free  In  the  sunlit  portion  of  Its  orbit. 
This  would  correspond  to  an  assembly  In  sunlight. 
At  time  t-0,  the  structure  moves  into  shadow  and 
deforms  until  reaching  a  new  thermal  equilibrium. 
Because  the  structure  is  In  a  geosychronous 
orbit,  the  time  to  cross  the  penumbra  is 
negligible  and  the  transition  Into  shadow  is 
considered  Instantaneous.  See  table  2. 
Furthermore,  all  Initial  cross-sectional 
temperature  gradients  are  neglected  and  the 
Initial  temperature  of  each  member  Is  Its 
radiation  equilibrium  temperature,  found  by 


3 


Discussion  of  Results 


For  both  load  cases,  examples  of  thermally 
induced  loads  versus  time  are  presented.  These 
curves  represent  the  Induced  mechanical  loads 
applied  throughout  the  structure.  Also,  the 
resulting  axial  stress  history  In  a  selected 
member  Is  given  as  a  ratio  of  maximum  axial 
stress  over  yield  stress.  This  stress  Is  further 
.'il  Into  bending  stress  and  stress  due  to 
extension.  For  Graphite  AS/Epoxy,  yield  stress  Is 
approximately  1  GPa  or  ISO  ksl. 


Case  1 

In  Fig.  7  the  axial  force  vs.  time  Induced 
In  all  members  using  cross-section  1  Is  given. 
As  expected,  the  magnitude  of  the  compressive 
force  Increases  smoothly  as  the  cross-section 
cools.  Because  a  uniform  Initial  temperature  was 
assumed,  the  Induced  moment  in  this  cross-section 
is  negligible.  Forces  and  moments  induced  In  the 
other  three  cross-sections  behave  similarly. 

The  resulting  axial  stress  for  such  a 
loading  Is  given  In  Fig.  8  for  structural  element 
number  4.  This  element  Is  located  at  the 
cantilevered  end  of  the  boom  and  Is  expected  to 
carry  a  higher  level  of  stress.  It  Is  evident 
that  stress  levels  are  not  overly  large  for  the 
structure  chosen.  At  a  time  corresponding  to 
that  required  to  cross  the  umbra,  the  maximum 
axial  stress  in  element  4  Is  20$  of  yield.  In 
addition,  stress  Increases  smoothly  with  the 
applied  load  Indicating  no  oscillatory  motion. 


Figure  7.  Induced  Axial  Force  In  Cross-Section 
for  Case  1 


Figure  8.  Ratio  of  Maximum  Axial  Stress  to  Yield 
Stress  in  Element  4  for  Case  1 


Case  2 

For  case  2,  both  axial  forces  and  bending 
moments  are  --oduced  by  the  thermal  loading.  For 
cross-section  1 ,  the  Induced  loads  are  given  In 
Figs.  9  and  10.  The  axial  force  Increases 
smoothly  and  approaches  Its  steady-state  value  at 
around  8000  seconds.  This  Is  then  the  time 


time  (seeiooo) 

Figure  9.  Induced  Axial  Force  In 
Cross-Section  1 
for  Case  2 

required  for  the  cross-section  to  regain  thermal 
equilibrium.  The  Induced  thermal  moment  behaves 
quite  differently.  An  early  peak  Is  reached 
within  the  first  500  seconds  due  to  the  lag  In 
the  heat  transfer  around  and  through  the  cross- 
section.  At  this  point  the  temperature  gradient 
Is  at  Its  largest.  After  a  peak  value  Is 
reached,  the  moment  decreases  to  its  steady-state 


value  at  8000  seconds.  The  results  from  the 
other  two  selected  cross- sect  Ions  differ  only  In 
magnitude. 

figure  11  gives  the  maximum  axial  stress 
induced  in  structural  element  4  by  the  loads 
above.  Steady-state  stress  values  are  on  the 
order  of  30*  of  yield.  Also,  the  stress 
increases  smoothly,  indicating  no  oscillation. 
It  is  interesting  to  note  that  the  stress,  and 
therefore  deformations,  are  in  phase  with  the 
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Figure  10.  Induced  Bending  Moment  in 
Cross-Section  1  for  Case  2 
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Figure  11.  Ratio  of  Maximum  Axial  Stress 
to  Yield  Stress  in  Cross-Section  1 
for  Case  2 

axial  forces  only.  This  indicates  that  axial 
forces  are  the  predominant  load.  One  then  might 
conclude  that  thermally  induced  bending  is  not 
significant  for  the  structure  and  load  chosen. 
Whether,  this  is  a  trait  in  most  structures  or  for 
different  load  cases  is  unclear. 


Conclusions 

This  study  has  attempted  to  investigate  the 
significance  of  thermal  loading  on  large  com¬ 
posite  space  structures  by  examining 
thermoelastic  responses.  For  this  purpose,  an 
Integrated  finite  element  model  for  performing  an 
uncoupled  thermoelastic  analysis  was  outlined. 
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temperature  gradients  within  members  and  for 
nonlinear  effects  such  as  radiation  boundary 
conditions  and  temperature  dependent  material 
properties.  Initial  studies  using  this  model 
were  carried  out  to  determine  the  significance  of 
thermally  induced  loads.  Although  only  limited 
studies  nave  been  carried  out,  several 
conclusions  concerning  thermal  loads  have  been 
reached. 

First,  the  present  research  indicates  that 
for  the  structure  modeled  herein,  there  is  no 
significant  dynamic  response  due  to  the  thermal 
loading  associated  with  entering  and  exiting 
earth’s  shadow.  Although  the  thermal  environment 
changes  instantaneously,  thermally  induced  axial 
forces  and  bending  moments  are  produced  In  a  much 
slower  ramp  type  fashion.  These  loads  require 
several  thousand  seconds  to  reach  their  steady- 
state  values.  It  has  been  estimated  that  even 
very  large  flexible  structures  will  have  their 
first  fundamental  frequencies  In  the  range  of 
0.01  to  10.0  hz.  The  fundamental  periods 
corresponding  to  these  frequencies  prove  much  too 
small  to  allow  excitation  from  the  Induced 
thermal  loads. 

Furthermore,  in  both  cases  studied,  the 
axial  forces  are  predominant.  However,  the 
thermally  induced  loads  are  very  dependent  on  the 
initial  reference  state  chosen.  Choosing  a 
reference  temperature  near  the  steady-state 
temperature  reduces  the  induce  axial  rorces  while 
leaving  the  bending  moments.  unchanged. 
Therefore,  it  is  not  yet  clear  if  thermally 
Induced  moments,  and  likewise  cross-sectional 
temperature  gradients,  can  be  neglected. 

Finally,  the  magnitude  of  the  thermally 
induced  axial  stress  is  less  then  30*  of  yield 
stress.  Although  this  cannot  be  neglected,  it  Is 
not  excessive.  Also,  these  stress  levels  could 
be  very  important  in  conjunction  with  stresses 
produced  by  manueverlng  and  docking. 

Although  the  study  presented  here  is  only  a 
preliminary  one,  several  important  trends  are 
evident.  However,  the  evaluation  of  thermally 
induced  axial  and  bending  forces  depends  on  many 
factors  and  will  require  more  investigation  in 
order  to  determine  their  real  significance  In 
thermoelastic  responses. 
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Abstract 


The  transverse  vibration  solution  for  both  uniform  and 
nonuniform  (variable  cross-section)  beams,  with  homogeneous  or 
nonhoraogeneous  (such  as  composites)  linear  elastic  material 
properties  has  been  widely  discussed  and  is  well  documented. 
However,  spatial  variations  in  material  properties  due  to  history 
dependent  damage  have  not  been  previously  considered.  In  this 
paper,  the  internal  state  variable  concept  is  Introduced  to 
account  for  the  cumulative  stress-induced  damage  in  a  simple 
composite  beam  structure  resulting  from  a  fatigue  type  loading  of 
several  hundred  thousand  cycles.  Furthermore,  the  transverse 
vibration  solution  based  on  the  history  dependent  material  is 
presented  for  a  simple  beam. 

Introduct ion 

A  number  of  analytical  approaches  to  the  well  known  partial 
differential  equation  for  the  flexural  vibration  of  a  nonuniform 
beam  structure  have  been  developed  since  Kirchhoff's  work  in 
1879[1].  He  solved  the  vibration  problem  for  wedged  and  cone- 
shaped  beams  in  terms  of  Bessel  functions.  With  the  invention  of 
the  digital  computer,  the  sophisticated  mathematical  approach  is 
replaced  by  numerical  approximation  without  loss  of  engineeriing 
accuracy.  Since  the  early  1960's,  numerous  approximate  solutions 
for  several  profiles  of  nonuniform  beams  have  been  documented 
using  the  finite  element  method[2,3]«  Most  of  these  solutions 


were  developed  for  beams  with  homogeneous  material  properties. 
Some  were  for  nonhomogeneous  (composite)  beams.  They  are  all 
restricted  by  the  assumption  that  the  section  properties  and 
stiffness  are  constant  in  time.  Neither  material  damage  nor 
environmentally  induced  degradation  is  considered.  The  use  of 
composite  materials  in  space  structures  makes  evaluation  of 
cumulative  damage  an  important  extension  for  beam  problems.  After 
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several  months  or  years  of  service,  it  is  expected  that  spatially 
variable  changes  in  material  properties  caused  by  cumulative 
load-induced  damage  will  create  new  complications  for  control 
engineers . 

Due  to  the  occurrence  of  load  induced  and  history  dependent 
damage  in  composites,  these  previously  obtained  solutions 
represent  inaccurate  approximations  of  the  actual  structure  after 
a  period  of  service.  In  particular,  the  resonant  frequencies  and 
mode  shapes  of  the  structure  can  be  severely  altered  by  the 
introduction  of  spatially  varying  damage  induced  stiffness  loss 
[4].  These  parameters  in  turn  will  have  a  substantial  impact  on 
the  active  control  algorithm  employed  for  control  of  flexible 
body  modes.  By  introducing  the  damage  parameter,  the  section 
properties  and  modulus  of  the  beam  will  vary  In  both  space  and 
time  due  to  the  history  of  stress  distribution.  The  stiffness 
change  can  significantly  shift  the  natural  frequencies  and  mode 
shapes.  With  the  damage  effects  involved,  the  differential 


equation  thus  becomes  difficult  if  not  impossible  to  solve  in 


closed  form. 

Experimental  research  on  advanced  composite  materials 
indicates  that  the  time  scale  for  damage  is  very  long  compared  to 
the  first  fundamental  frequency  of  the  structureC 5] .  Therefore, 
the  mathematical  algorithm  for  the  governing  equation  is  treated 
as  linear  with  slowly  varying  coefficients.  The  beam  equation 
thus  becomes  quasi-1 inear ,  which  can  be  solved  by  a  time 
integration  technique  at  each  time  incremant.  It  will  then  be 
possible  to  predict  actual  structural  response  for  a  fatigue 
type  loading  of  several  hundred  thousand  cycles  through  actual 
time  integration. 

In  this  paper,  particular  interest  is  being  placed  on  the 
natural  vibration  solution  of  a  simple  beam  structure  with 
spatially  and  time  varying  stiffness  caused  by  fatigue  type 
loading,  and  the  investigation  of  the  possible  effect  of  material 
damage  on  the  natural  frequencies  and  mode  shapes  of  the  beam 
with  various  boundary  conditions  (simply  supported, 
clamped-f ree ) .  The  concept  of  an  internal  state  variable  (ISV) 
is  Introduced  to  describe  the  history  dependent  damage  growth 
[8],  A  single  scalar  parameter  D  is  used  as  a  local  ISV 
representing  the  damage  state.  Together  with  the  internal  state 
variable  growth  law,  the  finite  element  solution  technique  is 
modified  to  account  for  the  time  dependent  stiffness  of  the  beam 


element . 


Problem  Summary 

In  order  to  Investigate  the  possible  effect  on  dynamic 
response  caused  by  material  damage  in  beam  structures,  a 
representive  heterogeneous  planar  beam  with  a  uniform  cross 
section  is  selected  as  shown  in  Fig.1.  With  the  applied  load 
history  along  the  transverse  direction  as  shown  in  Fig. 2  at  one 
of  the  nodes,  under  the  assumption  of  no  material  damping,  the 
beam  eventually  will  respond  in  the  first  fundamental  mode  which 
will  create  fatigue  type  cyclic  stress  distribution  for  the 
structure.  Finite  element  modeling  is  implemented  to  trace  the 
dynamic  stress  distribution  at  each  grid  point.  The  history  of 
peak  stresses  is  used  in  a  simplified  internal  state  variable 
growth  law,  which  thus  quantifies  the  spatially  varying  stiffness 
for  the  beam  element.  Using  a  spatially  varying  modulus  E,  a 
sensitivity  study  on  modal  analysis  is  conducted  to  show  the 
effect  of  material  damage  on  natural  frequencies,  and  mode  shapes 
of  the  heterogeneous  beam.  Details  of  the  finite  element  model, 
material  damage  model,  and  discussion  of  graphical  results  are 
presented  in  the  following  sections. 

Finite  Element  Model  with  Variable  Section  Properties 

Fig. 3  illustrates  the  simplified  heterogeneous  planar  beam 
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that  the  beam  is  composed  of  composite  material  with 


continuously  varying  section  properties  which  could  be 
discretized  as  shown  in  Fig. 3.  The  modulus  weighted  section 
properties  for  the  beam  at  any  time  are  defined  by[6] 
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where  E^.  and  A ^  are  the  modulus  and  cross-sectional  area  of  the 
ith  discrete  portion,  and  Eq  is  the  reference  modulus.  The 
effective  modulus  can  be  approximated  by 
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where  n  is  the  number  of  laminates  in  composite. 

Equatlon(4)  is  the  well  known  governing  equation  for  the 
transverse  vibration  of  an  undamped  planar  beam 

1x2( EI^x2) +m|t2-f (x,t)  (4) 


where  E  is  Young’s  modulus,  I  is  the  moment  of  inertia  of  the 
cross  section,  A  is  the  cross  section  area,  ra«  f A  is  the  mass  pe 
unit  length,  ^  is  the  material  density,  f  is  the  distributed 
force,  and  w  is  the  transverse  displacement.  If  we  consider  the 


spatially  varying  and  history  dependent  modulus  and  section 
properties,  the  governing  equation  can  be  rewritten  as 


2  12 
Z)  (J  W  fl  w 

flx2[E*(x,D,t)I*{x,D,t)^x2]+m(x)3t2-f(x,t)(5) 


where  the  parameter  D  is  an  internal  state  variable  which 
describes  the  damage  of  material,  and  t  is  time. 

As  mentioned  before,  the  scale  of  time  of  material  damage 
is  very  long  compared  to  the  fundamental  vibration  frequency. 
Thus,  the  governing  equation  is  treated  as  linear  with  slowly 
varying  modulus  and  section  properties.  Using  the  serai di screte 
variational  formulation  at  a  given  time  step[7],  the  equations  of 
motion  for  a  structural  element  become 


(6) 


where  and  are  the  element  stiffness  and  mass  matrix, 

respectively,  and  w  is  the  transverse  nodal  displacement  vector. 
That  is, 


E*(x,D,t)I*(x,D,t)N.N.dx 
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where  ^  are  the  shape  functions.  The  beam  is  modeled  by  four 
cubic  elements.  In  order  to  account  for  spatial  variation  in 
properties,  five  integration  points  are  used  to  evaluate  the 
stiffness  and  mass  matrices  of  each  element.  They  are  combined  to 
represent  the  coefficient  matrices  of  the  global  equations  of 
motion  for  transverse  vibration  of  the  heterogeneous  beam. 

[M]{w}+[K]{w}-{F}  (9) 

where  { F }  is  the  nodal  force  vector. 

The  Newmark  integration  technique  is  applied  to  solve  for  the 
dynamic  displacement  response.  After  solving  for  local  node 
forces  for  each  element  by  the  displacement  function  approach, 
the  modulus  weighted  section  property  is  again  used  to  calculate 
the  local  internal  stress  distribution  at  each  time  step: 

CJ^X  (x,y,D,t)*M(x,D,t)y/I*  (10) 

where  M  is  the  effective  nodal  moment,  and  0^*  is  the  axial 
stress  due  to  bending  .  Assuming  the  damage  occurs  at  the  peak 
stress  on  each  cycle,  the  stiffness  and  mass  matrix  are  updated 
whenever  the  grid  point  stress  reaches  the  peak  value.  Natural 
frequencies  and  mode  shapes  are  collected  at  a  period  of  every 
two  hundred  cycles  by  the  Jacobian  iteration  method. 


Material  Damage  Model 

Failure  of  continuous  fiber  composite  structural  components 
is  preceded  by  a  sequence  of  mi crostructural  and  macrostructral 
events  such  as  raicrovoid  growth,  fiber-matrix  debonding,  edge 
delamination,  matrix  cracking  and  fiber  failure,  etc,  which  are 
loosely  termed  damage.  Experimental  evidence  suggests  that  matrix 
cracking  will  predominate  prior  to  the  development  of  a 
characteristic  damage  state[8].  A  single  matrix  cracking  damage 
model  in  this  beam  problem  will  considerably  simplify  the 
solution  procedure  without  loss  of  engineering  accuracy. 

The  concept  o^  an  internal  state  variable  is  introduced 
through  constitutive  equations  to  quantify  the  damage  state  of 
matrix  cracking  density.  Due  to  the  history  dependent  reduction 
in  stiffness,  the  stress-strain  relation  is  of  the  form 
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..  is  thermal  strain,  C.,.,  is  the  effective  modulus 
■kl  ijkl 


tensor  which  is  given  by[8],  and  ,  (T^  are  the  second  order 


strain  and  stress  tensor,  r es pe c t i vely , 


f  ?  P 

C  i  j  kl  “  Cijkl" 


(12) 


where  ^3^  ^  is  a  set  of  material  constants, 


,  ^are  a  set  of 
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internal  state  variables  with  p  ranging  from  one  to  the  number  of 
damage  modes,  which  are  governed  by  the  following  ISV  growth  law: 

<*[  =  rL  (ifl  ,T,  0/t}  )  (13) 

In  the  one-dimensional  approximation  to  the  planar  beam 
problem,  a  single  s calar- val ued  internal  state  variable  is  used 
to  quantify  the  matrix  crack  density  of  the  composite  beam.  The 
constitutive  equation  thus  reduces  to 
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where  ,  £„  are  the  uniaxial  stress  and  strain,  is  the 
uniaxial  thermal  strain,  and  E'  is  the  local  axial  modulus  of  the 
beam  .  E'  is  related  to  internal  state  variable  D  via  the 
following  relation 

E*-Eo-^D  (15) 

where  Eq  is  the  undamaged  modulus,  D  is  the  ISV  which 
quantifies  the  matrix  cracking  density,  and  (3  is  a  material 
constant  which  is  determined  experimentally. 


Experimental  research  on  graphite  epoxy  composites  indicates  that 
the  matrix  crack  density  of  a  composite  specimen  under  constant 
fatigue  loading  can  be  represented  on  a  log-log  scale  by [ 9  3 


(16) 


Ln  D -K1 +  K2Ln  N 


where  N  is  the  number  of  cycles  of  fatigue  loading,  and  K1 ,  K 2 
are  material  constants  which  depend  on  the  stress  in  the  lamina. 
After  some  mathematical  manipulation,  an  internal  state  variable 
growth  law  is  derived 
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The  last  term  of  equation(17)  is  the  correction  term  for 
variable  cyclic  loading.  In  the  present  beam  problem,  the 
external  loading  is  held  constant  during  the  cyclic  response.  The 
interior  stress  will  decrease  due  to  the  loss  of  stiffness  caused 
by  material  damage.  Thus,  equation(l6)  represent  an  upper  bound 
of  the  damage  state. 

Equation(17)  completes  the  description  of  the  damage  model. 
Integration  of  this  equation  will  generate  the  current  value  of 
matrix  cracking  density,  which  is  substituted  into  equation(12) 
to  calculate  local  current  stiffness. 


Discussion  of  Computational  Results 


distribution  changes  due  to  the  difference  in  stress 
distributions.  For  each  case,  three(two  in  simply  supported) 
sample  points  were  chosen  to  trace  the  history  of  local  damage 
and  stiffness.  Natural  frequencies  and  mode  shapes  were  collected 
periodically  via  the  Jacobian  iteration  technique.  Assuming  pure 
bending  in  the  beam  structure,  no  damage  occurs  in  the 
neighborhood  of  the  neutral  axis.  Each  case  is  discussed 
separately . 

For  the  case  of  a  clamped-free  beam  the  maximum  stress  at 
each  sample  point  during  a  cycle  is  about  44.5  ksi  at  point  1, 
31.1  ksi  at  point  2,  18.5  ksi  at  point  3  (Fig.1).  Point  1 
undergoes  the  maximum  damage  after  one  hundred  thousand  cycles  of 
loading.  To  the  right  of  point  3,  no  significant  damage  occurs. 
Fig. 4  shows  the  local  damage  state  history  of  the  three  chosen 
points,  where  N  represents  the  number  of  cycles.  Fig. 5  shows  the 
stiffness  degradation  of  the  points.  Due  to  the  creation  of 
spatially  variable  damage  in  the  structure,  the  natural 
frequencies  decrease  as  shown  in  Fig. 6.  For  higher  vibration 
modes,  the  natural  frequencies  decrease  less  than  lower  modes. 
Although  the  maximum  stiffness  loss  is  at  point  1,  most  of  the 
beam  is  not  damaged.  Therefore,  there  is  no  significant  change  in 
the  first  five  fundamental  vibration  mod  -. 

For  the  case  of  a  simply  supported  beam  the  maximum  stress  at 
each  sample  point  during  a  cycle  is  about  44  ksi  at  point  1,  36.4 
ksi  at  point  2  (Fig.1).  Point  1  undergoes  the  maximum  damage 


r~~  * 


;  / 


Q  Q 
Q-  Q_ 


Fiq.  8  Local  Modulus  History  Cease  2) 


p  *  ■  i  'j 1  t  r*  ■?  wrmf 


rr 


i 


I 


K 


■  v»  vr»  \r"  i.  >>  -  \-v  ,-»  , 


after  one  hundred  thousand  cycles  of  loading.  No  significant 
damage  occurs  near  the  support.  Fig. 7  shows  the  local  damage 
state  history  of  the  two  chosen  points.  Fig. 8  shows  the  stiffness 
degradation  of  the  points.  Due  to  the  creation  of  spatiall 
variable  damage  in  the  structure,  the  natural  frequencies 
decrease  as  shown  in  Fig. 9.  Although  the  maximum  stiffness  loss 
is  at  point  1,  most  of  the  beam  is  not  damaged.  Therefore,  there 
is  no  significant  change  in  the  first  five  fundamental  vibration 
modes . 

Conclusions 

The  current  paper  has  attempted  to  develop  a  computational 
algorithm  for  performing  dynamic  analysis  of  composite  beams  with 
history  dependent  damage.  Because  the  algorithm  requires  that  the 
entire  load  history  be  considered  in  order  to  evaluate  the  damage 
state,  it  is  necessary  to  perform  the  numerical  scheme  over 
several  hundred  thousand  load  cycles.  This  procedure  is  very  time 
consuming,  requiring  approximately  22  CPU  hours  for  each  case  on 
a  Perkin-Elmer  3210  mini-computer.  Therefore,  the  implementation 
of  inelasticity  due  to  damage  presents  a  quite  cumbersome 
computational  task. 

The  results  for  the  two  cases  considered  herein  indicate  that 
although  the  fundamental  frequencies  of  the  structures  are 
significantly  affected  by  damage,  the  mode  shapes  are  not.  This 
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result  can  have  significance  with  respect  to  flexible  control 
algorithms.  For  example,  point  controllers  could  not  need  to  be 
replaced  at  some  other  point  on  the  structure  after  damage  has 
occurred.  On  the  other  hand,  previous  research[4]  has  indicated 
that  for  structures  such  as  trusses  which  are  primarily  loaded  in 
uniaxial  extension  and  compression,  the  mode  shapes  are 
significantly  altered  by  history  dependent  damage.  The  results 
would  seem  to  indicate  that  the  difference  in  those  results  and 
the  findings  here  are  that  when  bending  causes  through-thi ckness 
variations  in  damage,  the  effect  on  mode  shapes  is  significantly 
decreased.  It  should  be  pointed  out  that  while  the  model  used 
herein  assumes  that  bending  does  in  fact  produce  a  through- 
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